Modeling Cell Dynamics and Epigenetic Transitions. From Agent-Based to transport frameworks

Full text

Text extracted from the article PDF for reading and search. Tables, figures and equations are in the PDF, which remains the authoritative version.

1. Introduction

Cells are complex dynamic systems in which genetic information, environmental cues, and stochastic fluctu-ations interact in a nonlinear manner, controlling pro-cesses like proliferation, migration, apoptotic or necrotic death, cell differentiation, or tissue production/resorption [1]. Epigenetics refers to the varied set of molecular mech-anisms that regulate gene transcription without altering the underlying DNA sequence [2]. This regulation oc-curs through processes like DNA methylation, histone acetylation and methylation, nucleosome remodeling, chro-matin compaction, and large-scale reorganization of three-dimensional genome architecture [3, 4]. Together with DNA programming, they shape cellular identity (phenotype) and specific cell history (epigenetic markers accumulation), de-termining how a particular cell responds to external stimuli and stores information about past exposures [5, 6]. These mechanisms endow cells with a form of molecular mem-ory, which is partially reversible and heritable [2], thus contributing to cellular plasticity, a property central to many biological processes like embryonic development [7], immune adaptation [8], or cancer progression [9]. Predicting how these modifications evolve in time and space is therefore essential to understand how cells main-tain healthy tissue homeostasis as well as the progression and treatment of disease [10]. Given the complexity of these processes, with multiple spatial and temporal scales involved, computational modeling has become key to inter-pret and predict single cell and cell population evolution under varying environmental stimuli or therapeutic inter-ventions. Computational models allow virtual experimentation, connecting microscopic biochemical events to macroscopic observables such as tissue growth, shape conformation, or cellular heterogeneity. In particular, when referring to epi-genetic evolution, mathematical models help to elucidate how local enzymatic reactions and feedback loops in chro-matin can lead to global patterns of cell-state transitions [11, 12], or how environmental factors such as hypoxia and drug exposure reshape the epigenetic landscape of tu-mors, including cell heterogeneity and resistance acquisition [13, 14, 15, 16] Computational representations of cellular evolution typ-ically follow two paradigms: discrete and continuous. Dis-crete models represent individual cells as autonomous en-tities whose behavior is governed by probabilistic or rule-based dynamics [17, 18]. They naturally capture stochastic-ity, heterogeneity, and local interactions, and incorporate explicitly intracellular mechanisms such as gene regulatory networks or epigenetic circuits. However, they are compu-tationally expensive for large populations and often difficult to analyze mathematically. In contrast, continuous models describe the behavior of a whole cell population in terms of averaged fields, such as density, velocity, or mean internal state, that evolve in space and time under conservation laws and phenomenological constitutive relations [19, 20]. While these formulations lack internal molecular details, they are well suited for large-scale simulations, allow for systematic mathematical analysis, and link directly to con-tinuum theories of tissue biophysics thus permitting to model macroscopic systems.

A similar duality arises when modeling cell state evo-lution. Some formulations consider discrete transitions between a limited set of phenotypic states through stochas-tic or deterministic rules [21, 22]. While intuitive and easily tractable, such representations cannot reproduce the real gradual accumulation and erasure of epigenetic information [23, 24]. More realistic models introduce continuous epige-netic variables reflecting the smooth progressive transition of transcriptional shifts. This continuity gives rise to a land-scape of cell states that are not discrete but arranged along a continuous manifold in which cells occupy intermediate, metastable positions.

Within this continuous paradigm, two major strategies have been proposed. The first attaches to each cell or cell cluster in a representative volume (cell RVE) a set of internal variables that describe the epigenetic state and whose evolution depends upon the local environmental cues and kinetic laws, as well as the epigenetic state itself [13, 14]. The second considers the epigenetic variables as continuous distributions. In these models, each epigenetic variable is assimilated to a new independent coordinate, additional to the spatial ones, thus extending the physical space that now includes the epigenetic continuous manifold. The cell population evolves in this augmented spatio-epigenetic domain through transport and reaction equations [25, 15]. Although conceptually distinct, these two formulations describe the same phenomenon from different perspectives. In the independent-variable picture, the focus lies on the evolution of the population density and fluxes observed at fixed points in the extended space, highlighting the spatio-temporal evolution of the epigenetic state field. In the dependent-variable picture, one follows individual trajecto-ries of material points (cell RVE), carrying their internal variables as they evolve. Recognizing this equivalence al-lows for a unified interpretation of models that have been often considered as different in biological contexts. The main objectives of this work are therefore to review the main formulations of cell dynamics models from both continuum and discrete perspectives and to examine how phenotypic evolution is incorporated into these frameworks. Furthermore, we evaluate its conceptual and computational characteristics, establish the relations between the variables used in each approach, and identify the biological scenarios where each of them is more appropriate.

The paper is organized as follows. Section 2 discusses the discrete and continuous computational perspectives for modeling cellular processes without considering epigenetics. This is useful to introduce concepts and formulations that are extended in section 3 to discrete and continuous repre-sentations of epigenetic and phenotype dynamics. We also review different families of epigenetic models that can be found in literature. Under this light, we discuss in section 4 their respective advantages, disadvantages and limita-tions for each specific application. In this same section, we summarize the relations between the different variables that appear in discrete and continuous models and between Eulerian and Lagrangian approaches for the latter. Finally, Section 5 discusses the main conclusions and conceptual contributions of this work, and some perspectives for future research.

2. Modeling Cellular Processes Without Cell State

Changes Cells are living systems whose behavior emerges from the interplay of molecular reactions, gene regulation, bio-physical signaling pathways, and physical interactions with neighboring cells and the extracellular matrix. Processes like proliferation and death, differentiation, migration, shape change, or new tissue production/resorption are fundamental events that drive tissue homeostasis, growth, morphogenesis, and disease. Understanding and predicting how these processes evolve is therefore key for biological processes like regenerative medicine, cancer therapy, devel-opmental biology, or tissue engineering [26, 27, 28, 29] Two broad paradigms dominate physically-based compu-tational models of cell population behavior: discrete (agent-based) and continuous (field-based) approaches. Each embodies a different balance between mechanistic detail, computational cost, and analytical tractability. Finally, machine learning (ML) approaches complement the above by learning predictive mappings between inputs and out-puts directly from data [30, 31] or combining physics and data in hybrid frameworks such as Physics-Informed Neural Networks (PINNs) [32, 33]. Although the integration of purely data-driven models with classical biologically-based approaches promises to accelerate discovery and enhance predictive power in systems biology, in the following we shall focus only on the latter.

2.1. Modeling Cell Populations with A Numerable

Number of Cells Discrete models represent cells as individual entities that carry properties like position, velocity, or phenotypic state that evolve over time according to specific rules, with the aim of defining the value of each of those properties at each time instant. They also account for cell appearance and disappearance due to division or death, as well as entry into or exit from the domain.

Discrete cell-based models encompass a broad class of approaches, including cellular automata, lattice-based mod-els, and off-lattice individual-based formulations. Among these, we focus on Agent-Based Models (ABMs) [34, 35], which describe the coupled behavior of a set of agents, here identified with cells, that make decisions autonomously ac-cording to prescribed rules based on both external signals and their internal state. ABMs are suitable for modeling tumor evolution incorporating heterogeneity [36, 37] as well as epigenetic regulation and adaptation, as we will see in the next section. In particular, ABMs naturally capture cell-to-cell heterogeneity, allow for the incorporation of internal variables representing epigenetic states and cell state, and provide a flexible framework to model history-dependent behavior and memory effects at the single-cell level. Within this framework, each cell k is characterized by its position xk(t) ∈ Ωx ⊂ R3 that evolves over the simulated time interval I. Cell dynamics arise from key biological processes such as proliferation, death, migration, and phenotypic or epigenetic transitions. Cell proliferation and death are typically modeled through transition probabilities between discrete states (represent-ing, for example, different phases in the cell cycle), often described as Markov processes. Cell migration may be mod-eled either on a lattice, considering only a discrete number of sites that can be occupied according to probabilistic rules [38], or in a continuous space (off-lattice models) [17, 39], where cell movement is governed by a velocity vk that depends on both internal state and environmental cues. A common representation combines directed and random components:

where ζ is a unit vector with random direction, and u a unit vector indicating the direction of biased migration induced by the combined action of different taxes and drag forces, while A, B define the amplitude of these two components of the velocity, that may depend on the value of such taxes and the cell state. This equation must be solved simulta-neously for all cells, having therefore a discrete number of Ordinary Differential Equations (ODEs), associated withthe evolution of each cell k in space, which provides its location xk(t) ∈ R3 at time t.

Overall, cell evolution in ABMs (and discrete models in general) adopts a material description in which variables are associated to individual cells, enabling direct access to single-cell dynamics at every time t.

2.2. Modeling Cell Populations with A Continuous Cell

Den-Sity Continuous models take the complementary viewpoint by describing cell populations as continuous fields of av-eraged variables. These approaches, by their very nature, smooth out discrete events and individual fluctuations, so they are good to elucidate mechanisms of pattern forma-tion, tissue growth, and invasion waves in tumor invasion [28, 20], but fail when cell number is small or cell response is mainly driven by stochastic effects. Indeed, macroscopic continuous fields can be obtained from ABM models using mean field thermodynamics.

Continuous models are typically expressed in terms of systems of partial differential equations derived from con-servation laws. Two alternative but equivalent perspectives are used to describe the evolution of continuous systems: the Eulerian and Lagrangian formulations [40, 41]. Their conceptual distinction lies in the respective frame of ref-erence: while the Eulerian framework describes the time evolution of field variables in a fixed space, the Lagrangian framework follows individual small volumes of material (points in the continuum formulation) as they move and change over time.

Eulerian description. Using the general transport theorem, we can formulate the Eulerian differential version of the change rate for cell concentration, as (2)

where D is the so-called material derivative or total deriva-tive with respect to time, that omitting the dependencies for the sake of clarity is expressed as (3)

and v(x, t) the spatial advection velocity of cells at point x and time t that depends on the external signals (taxes) and the state of the cells, becoming clear that we need to assumethat all cells in a cell RVE (x) have the same velocity. Also, rp is the rate of proliferation (rate of new cells production per unit cell) at x and t, and subjected to environmental signals µ(x, t) (whose evolution may also be described with a reaction-diffusion equation), while rd is the rate of death by necrosis or apoptosis. These values can be modeled with different degree of detail, from just predefined functions fitted from available data, to very complicated models encompassing the intrinsic cell or molecular mechanisms involved [42]. The term qd is the diffusive flux of cells, that is, cells entering the reference volume through its surface (cells per unit time and unit surface) due to diffusion or pedesis (random movement of cells) and is an emergent term appearing from thermodynamic averaging. Thus, the explicit dependencies of each variable that we will omit from now on for clarity purposes are the following:

We can express the advection flux as qa = cv, and using divx(cv) = c divx(v) + v · grad c, Equation (2) can be rewritten in the conservative form as:

where the source term is r = rp − rd.

If we assume a simple linear diffusion model, we have, in the end, (5)

being Dc the diffusion tensor.

We therefore need to make explicit the diffusion matrix Dc at each point x as well as the reaction function r, and the advection flux qa. This latter may be obtained after defining or computing two velocities. First, the advection velocity has to be determined at each point (e.g. solving a fluid flow or solid deformation problem to get the dis-tribution of drag forces or velocities acting on the cells at every spatial point and time due to the surrounding solid or fluid velocity). Second, a phenomenological expression is required for the intrinsic cell movement due to external signals (taxes). If we also want to model the evolution of the concentration of these signals, such as oxygen or other nutrients, this can be expressed in terms of transport equations for each species or signal µi, which in Eulerian form write:

where, here, a simple linear diffusion is considered. Finally, si is the source term for species i that comprises decay as well as production and/or consumption by cells, being expressed as:

with di the decay rate, and αi the function modeling the net production/consumption rate of each species by cells, located at point x at time t.

This context is well suited for describing flows, diffusion, and reactions making abstraction of the actual evolution of properties of the individual cells composing the system. It is, therefore, particularly useful for problems where the precise location of individual cells is less relevant than the spatiotemporal evolution of key concentrations. Such an approach is ideal for studying the emergence of biological patterns, as seen in morphogenesis [43], tumor growth [20], and the dynamics of heterogeneous cell populations [44, 45]. However, we cannot follow directly the history of the variables associated to each cell (velocity, position, state,..) since the cells passing through a certain fixed point x change over time. This makes it more complex to account for history-dependent behaviors, particularly in processes driven by the evolution of the epigenetic state. While certain techniques can mitigate this by incorporating additional state variables, they significantly increase model complexity; therefore, a Lagrangian description is often more suitable for such cases. Lagrangian description. The Lagrangian approach focuses on following each material point, usually defined by its coordinates X at an initial time t = t0, along its trajectory through space and time. Rather than examining what hap pens at fixed spatial coordinates, one tracks how properties attached to each cell RVE, such as position, velocity, cell concentration or internal state, change as it moves and evolves. In this case, the initial position of each point of the domain of interest must be explicitly defined, like in the discrete case. To derive the governing equations in this framework, we first establish the motion mapping x = χ(X,t), which links the reference configuration to the current position and is obtained from the definition of the velocity field for all material points V (X,t) = v(χ(X,t),t). The deformation of the domain is characterized by the deformation gradi ent F = ∂x ∂X and its determinant, the Jacobian J = detF, which represents the local volume change. The material cell concentration is then defined as C(X,t) = Jc(x,t), repre senting the number of cells per unit of material reference volume. By applying the Reynolds Transport Theorem in its differential form, the Eulerian conservation law (Eq. 4) is pulled back to the reference configuration. In this de scription, the material derivative simplifies to the partial derivative with respect to time, D Dt = ∂ ∂t, and the advective terms are implicitly handled by the motion of the coordi nate system. The key identity, which follows directly from the definition C = Jc and the material time derivative of J:

allows us to identify the left-hand side of (4), with qa = cv, after separating the advective and diffusive fluxes, as:

Multiplying (4) by J and applying the Piola identity J divx b = DIVX(JF−1b) to the diffusive flux:

where the material diffusive flux is defined as Qd = JF−1qd. The advective terms are implicitly handled by the motion of the coordinate system, and the balance equation becomes:

where the material counterparts of c,rp,rd and qd, defined per unit volume of the initial configuration are:

The material flux Qd accounts for the transformation of the diffusive flux vector from the deformed spatial configuration back to the reference state. If we assume a simple linear diffusion model, using the transformation of the spatial gradient, gradx c = F−TGRADX(C/J), and substituting into the definition of Qd:

which can be written compactly as:

the material diffusion tensor. Substituting into (8) yields the closed material balance equation:

where R = Rp−Rd is the net material reaction rate. Analogously, we can derive the Lagrangian equation for the evolution of the chemical species. We define the material concentration of the i-th species as Mi(X,t) = µi(χ(X,t),t). Unlike the cell concentration C = Jc, the species concentration µi is an intensive quantity, so its material counterpart carries no Jacobian factor. Expanding the divergence in (6) using the product rule:

The first two terms on the left-hand side are the material derivative of µi, which equals ∂Mi/∂t by definition of the pull-back. Using divx v = ˙J/J, the equation becomes:

Applying the Piola identity to the diffusion term and using gradx µi = F−T GRADXMi:

where DM,i = F−1DiF−T is the effective material diffu sion tensor for the species Mi, consistent with the definition provided for the cellular phase. We thus obtain the conser vation law in the reference configuration:

where the material source term Si(X,t) = si(χ(X,t),t) accounts for the net production or consumption rate trans formed to the reference state. Substituting the Eulerian source term (7) and using c = C/J:

As in the Eulerian case, we need to define the reaction functions, Rp,Rd, Si, as explicit functions of the cell den sity C and the external signals µi. Since some of these functions depend on the value of the signals at the current spatial position x = χ(X,t), we also need to compute the velocity of each RVE, which depends on the drag forces and cell movement driven by external taxis signals already mentioned for the discrete models. By integrating this velocity:

we update the position x = χ(X,t), the deformation gra dient F, and the Jacobian J at each time step. The resulting equations in both Eulerian and Lagrangian cases are usually solved by numerical methods, such as Fi nite Elements or Finite Volumes, which use interpolations or approximations of the solution in terms of a number of discrete values in space. This discretization transforms the original partial differential equations into a system of ordinary differential equations (ODEs) with respect to time, a structure similar to that of the discrete case, al though with the flexibility provided by the huge amount of approximation methods available. While Eulerian methods typically employ a fixed spatial mesh, Lagrangian procedures offer more flexibility in track ing material deformation through two main approaches: the Total Lagrangian formulation, where all variables are referred to the initial configuration at t0 = 0; and the Updated Lagrangian procedure, where the reference con figuration is reset to the current state at the beginning of each time interval [t,t + ∆t]. In the latter, the reference coordinates and state variables are incrementally updated (Xt ←Xt+∆t) after each successful time step. As commented, the Lagrangian framework is widely used when one is interested in the history of individual elements [46, 47], and knows the initial position of each point. Consequently, it may be less suitable for problems where it is difficult to establish the initial configuration and/or the initial spatial distribution of each variable. 3. Modeling epigenetic evolution Cells change their phenotypic state by epigenetic changes induced by their microenvironment. Epigenetic regula tion shapes the long-term dynamics of cellular populations by enabling gradual and partially heritable modifications in gene expression without modifying the DNA sequence [48, 49]. Unlike genetic mutations, DNA methylation or histone modifications are reversible and can vary smoothly, generating a wide spectrum of cell behavior facilitating adaptation to environmental stresses, differentiate into spe cialized lineages, or acquire drug resistance [50, 51]. Over the last two decades, several families of compu tational models have been proposed to represent these processes. One first distinction is, again, between those that describe epigenetic dynamics as transitions between a discrete number of states and those that consider evolution in a continuous state space. 3.1. Modeling cell populations with a discrete number of phenotypes In these models, each cell can adopt a finite number N of distinct phenotypes {ϕ1,ϕ2,..,ϕN}, such as differentiated versus stem-like states or active versus repressed chromatin configurations. The dynamics of each cell phenotype can be formalized in any of the frameworks described in the previous section: namely ABMs or continuum equations in either a Lagrangian or Eulerian configuration. The ABM approach tracks the state Ek of each individual cell k. The probability of a state change during a time step ∆t is given by: (15)

This Markovian perspective is of course widely used as it allows for a tractable mathematical description of cell-state transitions and a direct interpretation of lineage-tracking data and single-cell plasticity. Alternatively, one can track the evolution of the popula tion density ci(x,t) for each phenotype ϕi at a fixed spatial location. The dynamics are then governed by a system of coupled flux balance equations following the ideas of Section 2:

where vi is the velocity field of phenotype i, Dc,i is the diffusion matrix and ri represents the net growth (prolif eration minus death). This Eulerian formulation can be mapped onto its Lagrangian counterpart via a pull-back op eration, resulting in a material balance equation consistent with (8). While the choice of framework determines how the col lective dynamics are computed, the underlying biology is governed by the transition rates or probabilities between states, which represent the cellular response to multiscale stimuli. It is important to emphasize that this kinematic choice—whether formulated in the current or reference configuration—is a matter of mathematical representation for the bulk motion and does not dictate the underlying phenotypic dynamics. The general transition rate κij from phenotype i to j (Eqs. 15, 16) allows us to categorize existing models based on which components are explicitly modeled and can be defined as a function of intrinsic stabil ity, environmental cues (denoted from now on collectively as µ in the spatial frame and M in the material frame), and cellular memory:

where Gij represents the intrinsic stability of the phe notypic attractor (determined for example by the Gene Regulatory Network (GRN) architecture), and Hij(x,t) is the integrated signal history (memory), defined as:

where Wij is a history matrix kernel. It must be noted that to compute signal history it is unavoidable to resort to a material description of the stimuli, in order to track the history of cell trajectories. Note that if we neglect history dependence, that is, Wij(u) = δ(u)I, then we obtain a pure Eulerian expression for κij:

Instead of computing the history term Hij, an equiv alent approach is to use internal variables. For a generic spatial vector ξ(x,t) defined as ξ = H(x,t), using Leibniz derivation rule (20)

It is important to note that as the cell movement depends on epigenetic changes, the previous equation is a nonlinear PDE on ξ illustrating explicitly why the Eulerian approach is not the natural one for epigenetic changes tracking. Now let us define Ξ(X,t) = ξ(χ(X,t),t). Then, (21)

This provides a material description of the internal epi genetic state. Next section explores this idea from the kinematic point of view. Evaluating Eq. (21) explicitly, the material evolution equation for Ξ reads:

The literature on phenotypic transitions can be broadly categorized into two families based on the primary driver of the switch. On the one hand, models focused on intrinsic dynamics and stochasticity, rooted in GRNs, describe tran sitions as a result of molecular feedback loops that create bistability or multi-stability [52, 53]. In these frameworks, transitions are often triggered by molecular noise, allowing cells to jump between stable attractors in an epigenetic landscape. Later extensions have focused on cancer and reprogramming [54, 55] or, at the single cell level, for mod eling the activity of chromatin regulators [56] using Markov processes. Alternatively, the population-level consequences of these transitions are often analyzed through Matrix Pop ulation Models (MPMs) [57, 58], focused on on asymptotic behavior (by analyzing the spectral properties of the stable phenotype distribution towards which the population con verges can be predicted) and applied to the study of how phenotypic plasticity (state) and other traits such as age af fect population dynamics in plants [59] or invertebrates [58]. This perspective is further developed in mesoscopic frame works where discrete traits are acquired through mutation probabilities [60], as well as in ODE models modeling pro cesses such as drug resistance with sensitive and resistant compartments [61] or the epithelial-mesenchymal transition (EMT) where fixed transition rates are fitted to genomic data [62]. In this class of models, focused on intrinsic dy namics and stochasticity and thus exploring the interplay between the internal network architecture and noise, the transition rate can be set as κij ≈ F(Gij) + ζ(t), where ζ(t) accounts for stochastic molecular noise. Transitions occur when stochastic fluctuations ζ(t) allow the system to overcome the potential barrier defined by Gij. Note that this random term has not been included in the general def inition of the transition rate in Eq. (17) to stay under the deterministic framework. Instead, the PDE framework cap tures the resulting mean-field transition rate, representing the collective probability of these stochastic jumps across the population. On the other hand, a second family of models explicitly accounts for external stimuli—such as mechanical, biochem ical, or social cues—as drivers of adaptation. In mechanobi ology, for instance, volumetric strain and fluid shear stresses are modeled as governing factors for cell differentiation dur ing bone fracture healing [63]. In mechanistic bone remod eling, mechanical strain regulates the local concentrations of signaling molecules such as RANKL, OPG, and TGF-β, which subsequently drive phenotypic changes [64, 65], while in wound healing, growth factors like PDGF and TGF-β act as the primary controllers of cellular state transitions. These external drivers often interact with cellular memory, as seen in hypoxia models where cells returning to normoxia do not revert to their original state but adopt persistent aggressive phenotypes [66, 67]. Another widespread ap plication of this approach is the study of drug resistance, typically employing two-compartment models for sensitive and resistant populations with a transition dependent on (22) drug exposure using ODEs [68, 69, 70] or stochastic models such as Markov branching processes that allow for random mutations [71, 72]. To bridge the gap toward continuous frameworks, these approaches could be extended by al lowing an arbitrary number of phenotypes arranged along fitness axes [73]. Furthermore, this environmental influence extends to evolutionary game theory, where phenotypic transitions are modeled as frequency-dependent strategies. In these cases, the fitness and transition probability of a specific phenotype are conditioned by the composition of the surrounding population, capturing the competitive or cooperative nature of cellular ecosystems [74, 75].

In this family, the transition is driven by external cues, leading to κij ≈ F(µ(x,t)). In these cases, the internal biology Gij is often treated as a constant threshold that µ must cross to trigger differentiation. In case of irreversible transitions [66, 67] the intrinsic term Gi acts as a one-way barrier with κij > 0 and κji = 0.

3.2. Modeling Cell Populations with Continuous

Phenotype A key challenge in modeling cellular adaptation and differentiation lies in representing the continuous evolu tion of internal state variables like chromatin organization, transcriptional activity, and metabolic configuration. In contrast to discrete models, continuous frameworks describe the epigenetic state of a cell by a set of real-valued variables. This approach aligns with the physical and biochemical nature of epigenetic modifications, which involve gradual addition and removal of biochemical markers. Two main families of methods have been recently pro posed within this framework: the first considers phenotype structured models in which phenotypic variables are mod eled as additional coordinates (independent variables), and the second in which phentotype is described by a set of dependent variables that evolve through local evolution laws without increasing the system’s dimensionality.

3.2.1. Local Continuous Distribution of The Epigenetic

Vari Ables: Structured Models Traditional continuum models of tissue dynamics op erate solely in the physical space, Ωx, describing how cell density or velocity evolves in time and position. However, now we want to consider that cells composing a certain RVE (or material point in continuous physics) may have different epigenetic states. If we want to track pointwise the evolution of this state distribution within a bounded set Ωξ, driven by epigenetic modifications or transcriptional re programming, we must extend the physical space to include additional dimensions that capture these new epigenetic coordinates. The state of the system is then characterized by the evolution of each property like cell density c(x,ξ,t), where x ∈ Ωx ⊂ R3 represents the physical position, and ξ ∈ Ωξ ⊂Rn represents the cell state [25, 76, 77]. The evolution c is governed by a transport–reaction–diffusion equation that expresses the conservation of the cell number in the spatio-epigenetic extended domain. Following (6), we can write again the Eulerian description as:

where we have removed the dependence on each variable since it is stated in (2). The reaction term rp−rd represents the net production of cells at the spatial point x, with state ξ and at time t, by proliferation or death. The flux qa ξ = cvξ encodes the advection movement of cells in the epigenetic space, com ing from transitions between epigenetic states induced by external signals. For instance, hypoxia or drug stress may induce a drift in ξ toward resistant or stem-like configura tions [15]. The diffusion term divξ Dξ gradξ c captures stochastic gene-regulatory fluctuations in the epigenetic state [21, 22], producing, for example, natural decay driv ing cells back towards baseline states in the absence of external signals. Boundary conditions at the edges of Ωξ limit the epigenetic transitions, effectively constraining the admissible phenotypic range. Therefore, the resulting cell distribution c(x,ξ,t) describes not only how cells move in the physical space, but also how their identity is distributed in the epigenetic landscape. Eq. (23) can be written using a compact notation as:

with vy = (vx,vξ)T and Dy = diag(Dx,Dξ). The dis tribution of cells in the real spatial domain ρ(x,t) can be recovered by:

To analyze history-dependent behavior and phenomena such as cell conditioning, we may transition from the Eule rian description presented above to its Lagrangian coun terpart. We define the material coordinate Y = (X,Ξ), representing the initial position and epigenetic state of a cell cohort at t = 0. The current configuration y = (x,ξ) is determined by the motion mapping y = χ(Y,t), which is obtained by integrating the coupled system of ordinary differential equations:

subject to the initial conditions x(X,0) = X and ξ(Ξ,0) = Ξ. Note that while the velocities vx and vξ depend on the current state (x,ξ) to reflect local microenvironmen tal stimuli (e.g., chemotactic gradients or stress-induced reprogramming), the resulting trajectories are uniquely identified by their material labels (X,Ξ).

In this framework, the material density C(Y,t) repre sents the number of cells per unit of initial spatio-epigenetic volume. It is related to the spatial density c through the Jacobian J = det(FY ), such that C(Y,t) = Jc(χ(Y,t),t). We can obtain the material analogous of (11) from Eq. (24). Applying the Piola identity in the extended space, J divy b = DIVY(JF−1 Y b), together with the pull-back of the extended spatial gradient, grady f = F−T Y GRADY f, and the block-diagonal structure of the spatial diffusion tensor Dy to obtain:

where R = Rp−Rd is the net growth rate evaluated along the trajectory and the effective material diffusion tensor in the extended space is:

It is important to note that even though Dy is block diagonal, DY is not in general, since the off-diagonal blocks of FY introduce cross-coupling terms between spatial and epigenetic directions. These blocks are:

which are generically non-zero: the first encodes how the current spatial position of a cell depends on its initial epigenetic state (cells with different phenotypes migrate differently), and the second encodes how the current epige netic state depends on the initial spatial position (cells in different microenvironments undergo different epigenetic evolution). Only when both blocks vanish does FY become block-diagonal, J factorizes as JXJΞ, and Eq. (26) reduces to (27)

A fundamental advantage of this Lagrangian approach is that the complex advection terms from the Eulerian de scription— namely divx(cvx) and divξ(cvξ)—are implicitly handled by the trajectory mapping χ. In the material frame, the coordinate system moves with the cell cohort and thus, the flow of cells through the epigenetic landscape is embedded in the evolution of y(t). Moreover, this for mulation effectively preserves the identity of cell lineages. This path-dependency is critical for modeling how prior environmental exposure determines the current adaptive capacity of tumor cells. The phenotype-structured framework has been exten sively applied to capture clonal selection and population heterogeneity in cancer evolution. Early foundational mod els in this area primarily focused on selection-driven dy namics, where phenotypic drift is neglected (vξ = 0). This is classical mean field theory obtained from random fluc tuations of epigenetic noise. In these cases, the evolution of c(x,ξ,t) is governed by differential proliferation and mi gration rates along the phenotypic dimensions and random mutations, modeled through diffusion in Ωξ. This approach has been applied to study the interplay between cells and oxygen levels [78, 16], the go or grow behavior of tumor cells [79] or drug resistance [76, 80]. More recently, the inclusion of an advective phenotypic drift (vξ ̸ = 0) has allowed to model microenvironment-driven plasticity, by adding a drift as in Fokker-Planck equation. Examples include [25, 15], where cells are driven to stem-like pheno types when exposed to hypoxia (see some results in Figure 1), or [81], where the Epithelial-to-Mesenchymal transition is modeled with phenotype changes driven by the density of the ECM. Additionally, in [82] the drift is caused by an epigenetic drug.

Fig 1: Results of cell evolution in the extended spatial phenotypic space using an structured model. The first and second rows respectively reproduce two cases and each column corre sponds to a time step of cell evolution along the phenotypic dimension s (where s = 0 corresponds to differentiated cells and s = 1 to cancer stem cells) and space X. Reproduced from [15].

Besides cancer, there are other applications of this ap proach ranging from general formulations under the prism of evolutionary biology for populations under fluctuating environments [81] to particular examples of bacterial persis tence under antibiotic stress, developmental differentiation or immune adaptation. The reader is referred to [83] for a more thorough review.

This extra-dimension formalism bridges the gap be tween microscopic stochastic models of molecular regulation and macroscopic tissue-scale descriptions, and is a versatile tool for representing continuous phenotypic plasticity in cell populations.

3.2.2. The Internal Variable Approach in Continuous

Epigenetic Evolution Inspired by the fundamentals of internal state theory in control theory [84] and continuum mechanics [85], an alter native approach to the ones described above was presented and applied in the work of Pérez-Aliacar and colleagues [13, 14], using internal variables to describe the phenotypic state of cells.

The internal variable formalism originates from non equilibrium thermodynamics, where it represents hidden degrees of freedom that store history-dependent informa tion [86, 87]. In the biological context, the internal variables represent the epigenetic state, encoding memory of past conditions such as hypoxia or therapy exposure. Therefore, we consider a finite dimensional spatial continuous distribu tion of epigenetic states, with values at each spatial point x and time t, as a real-valued finite-dimensional vector ξ. These internal variables are considered as properties that are carried by each cell RVE (each material point). They change along their trajectory under environmental stimuli, so their evolution is mathematically defined by an evolution equation similar to that of cellular concentra tion. Inheritance at division—full, partial, or noisy—can be represented by appropriate rules to conserve or dilute epigenetic memory across generations, providing a mini mal mathematical expression for epigenetic heritability in proliferating populations [14].

As the general framework formulated in [13], this ap proach considers m cell populations characterized by their densities ck(x,t). Each population k is endowed with dk specific epigenetic internal variables, ξk i(x,t) ∈ Γi,k ⊂ R (i = 1,..,dk), representing the epigenetic accumulation per unit cell within a RVE. Finally, the environment is described by n chemical species acting as signals, with mass concentrations µj(x,t) (j = 1,..,n). The evolution of these variables is derived from con servation laws [13]. The conservation equation for the epigenetic state per unit volume, ckξk i, in the absence of a diffusive flux for the epigenetic state (consistent with the assumption that epigenetic transport occurs solely through cell convection), reads:

where qk = ckv +qk d is the total cell flux (advective plus diffusive) and fk i is the epigenetic source term per unit cell. Expanding (28) using the product rule, substituting into (28) and collecting the terms multiplying ξk i :

The bracketed term is recognised as the left-hand side of the cell conservation equation (4):

where rk = rk p −rk d is the net proliferation rate. Substituting:

Solving for ∂ξk i/∂t and dividing through by ck yields the evolution equation for the epigenetic variable:

The convection term ensures consistency within the Eulerian framework as the state is transported with the cellular movement. Finally, fk i is the source term, expressed as:

where ωk i = pk i(µ(x,t),ξ(x,t)) models the acquisition of epigenetic changes, that might be induced by a combination of external stimuli, their rates and gradients, and the cur rent level of the internal variables themselves, which allows to model cell adaptation as a response to environmental changes. λk i is a decay coefficient representing spontaneous repair mechanisms. The third term accounts for phenotypic inheritance during cell division (governed by the growth rate rp); here, βk i represents the fraction of the phenotypic state inherited by daughter cells. In this regard, the term is proportional to the growth rate in cells rp through the coefficient βk i ∈ [0,1] that regulates inheritance in such a way that when βk i = 1 daughter cells inherit the phenotypic state of their progenitor, whereas when βk i < 1, only a per centage of daughter cells inherits the state of their parents. Finally, the fourth term incorporates phenotypic state de struction due to cell death and therefore it is proportional to death rate rd. By plugging Eq. (30) in Eq. (29), we obtain the final expression for the evolution of the internal variable (31)

The governing equations for cell populations and species are those detailed in the previous section, with the nuance that now the internal variables ξk = (ξk 1,..,ξk dk ) modu late the populations’ behavior by influencing its flux qk, growth and death rates, rk p and rk d, in the case of cells, and consumption/production αi in the case of species. That is, the internal variables representing the phenotypic state modulate the populations’ behavior thereby closing the coupling between the environment and cellular response. Finally, to accurately capture complex, biological pro cesses in which history-depencence is sought to be tracked individually, a Lagrangian formulation is often better suited than its Eulerian counterpart. In this framework, we track individual RVEs as material points, integrating the en vironmental stimuli they experience along their specific trajectories.

As in the previous section, let X denote the refer ence (initial) position of a material RVE at t = 0, and x =χ(X,t) its current position at time t, with the map ping both configurations obtained by integrating the mate rial velocity V (X,t), determined by the local microenvi ronmental signals and the current phenotypic state of the RVE.

Let Ξk i(X,t) = ξk i(χ(X,t),t) represent the epigenetic internal variables associated with the material point X per unit cell. By moving with the cellular RVE, the explicit convective term vanishes, and the conservation law (28) for the epigenetic state can be rewritten in its material form. First, we apply the pull-back of the three fields involved. The spatial gradient of ξk i transforms via the inverse transposed deformation gradient:

The cell concentration satisfies Ck = Jck, and the diffusive cell flux transforms according to the Nanson formula for vectors, Combining these three relations:

Expanding the scalar product and using FTF−T = I:

Substituting this result into the material form of (29), we obtain the exact Lagrangian evolution equation for the epigenetic variable:

where Fk i (X,t) = fk i (χ(X,t),t)−Rk(X,t)Ξk i(X,t) is the material source term. The advective part of the convective term, v·gradx ξk i, is absorbed into the material time deriva tive and thus disappears in the Lagrangian frame. However, the residual transport term 1 Ck Qk d · GRADX Ξk i persists: it accounts for the redistribution of the epigenetic state driven by the random (diffusive) motion of cells relative to the mean velocity field v, and vanishes only if this diffusive cell flux is negligible. Eq. (32) is analogous to Eq. (22). Note that, unlike the formulation for cellular concentration C, the Jacobian J = det(F) does not appear explicitly in the diffusion term for Ξk i. This is because Ξk i is defined as an intensive property per unit cell, making its transport independent of volumetric dilution. The term Fk i is the material source term, which tracks the net rate of phenotypic changes per unit cell as it navigates the domain:

where now Ωk i = Pk i (M(X,t)),Ξ(X,t))). In this setting, the internal variables act as a biological memory, accu mulating the effects of past exposure to stressors such as hypoxia or drug concentrations. This Lagrangian frame work bridges the gap between discrete cellular histories and continuum tissue mechanics, providing a robust pathway to track the emergence of drug-resistant niches in evolving glioblastomas.

In [13], this general framework (in its Eulerian for mulation) is particularized to the study of glioblastoma adaptation to hypoxia through several simplifications, i.e., there is only one cell type and one epigenetic variable, that evolves depending on hypoxia exposure without considering random mutations. As an illustration, figure 2 reproduces cell and phenotypic state evolution for a particular case. Furthermore, the internal variable approach has proved its ability to replicate in vitro data of drug resistance in glioblastoma spheroids [14] in this case using an ODE model with again one phenotypic variable, demonstrat ing that macroscopic observables (growth curves, viability, rebound after drug holidays) can be reproduced with a unique averaged value of the phenotypic variable at each point.

Fig 2: Cell and phenotypic state evolution with the in ternal variable approach. Example of GBM cells under an in silico cyclic hypoxia experiment, where cells migrate towards the oxygenated area while undergoing phenotypic changes, ranging from GBM differentiated cells to Cancer Stem Cells (CSCs) due to hypoxia exposure. Reproduced from [13].

To close this section, Table 1 summarizes the main features of the two families of models herein discussed for modeling phenotype as a continuum variable.

4. Discussion. Comparison and Conceptual Integration of

the Different Models The choice between these different paradigms (discrete versus continuous, structural versus internal variables, and Eulerian versus Lagrangian) depends on the data available, the scale, and the question of interest among other factors. In the following, we highlight the main characteristics of each approach versus its alternative, identify the main pros and cons of each of them for specific applications and establish the connections between them which facilitates an easy translation from one to another.

4.1. Discrete Versus Continuous Cell Models

From all the discussed in previous sections, it becomes clear that discrete models are attractive due to their simplic ity and biological interpretability. They provide a natural description of the system, with real stimuli and responses, catching nontrivial dynamic behaviors such as hysteresis, multistability, and oscillations [21]. They also allow direct correspondence with experimentally observed phenotypic categories, and their parameters—transition rates, pro liferation rates, or death rates—can be estimated from population-level data. Since they work at the individual scale, they allow to analyze how local interactions give rise to global behavior, thus capturing emergent phenomena at the population level [88]. Moreover, they inherently incorporate stochastic phenomena, allowing to control and analyze predictions uncertainty [89]. Finally, they have a modular structure and are much flexible, allowing to easily introduce new agents or rules, without modifying the entire model. However, they have a high computational cost due to their combinatoric nature when considering cell interac tions, which makes them prohibitive for large numbers of cells as happens when solving tissue-scale domains. Even more dramatic, they have many parameters, which in many cases translates into a lack of parameter identifiability [90] making the fitting task extremely complicated [91]. Finally, they tend to be more dependent on their implementation than continuous approaches [92].

Conversely, the continuum formalism provides direct access to distributions and fluxes that can be compared with experimental histograms obtained from flow cytome try, single-cell sequencing, or imaging. Furthermore, they naturally accommodate heterogeneous environments, and easily show how local environmental variables (e.g. oxygen, drugs, nutrients) influence on the drift and diffusion coeffi cients, leading to spatially dependent adaptation patterns [93, 76]. They also permit an easy derivation of macro scopic integral quantities—such as mean phenotype, total cell number, and spatial fluxes over surfaces, and, there fore, an easy coupling with continuum models of nutrient diffusion, mechanical stress, and therapy delivery. How ever, continuum formulations do not permit to distinguish between cells located at the same point with distinct prop erties if they are not considered as different entities, which implies an increase in the number of variables to model. New hybrid models are proposed with the aim of com bining the best of both worlds. For example, ABMs may be used to inform the parameters and closures of continuum equations [94, 95], while continuum models can provide averaged fields that guide stochastic simulations [96, 97].

Table 1: Comparative analysis of continuum phenotypic modeling frameworks: Phenotype-structured models and Internal variables. Feature Phenotype-Structured Internal Variables Mathematical Basis PDEs (Advection-Diffusion) in extended space (x,ξ)(x, \xi). PDEs in real space (x)(x) for cells and phenotype ξ\xi. Phenotype (ξ\xi) Continuous independent variable. Property attached to the point and connected with cells. Variables describing phenotype Infinite (continuous distribution). Finite (one value per internal variable). Computational Cost High: Requires mesh in n+1n+1 dimensions. Moderate: Standard mesh; extra DOF per node. Cell History Difficult to track in Eulerian grids (mixing). More natural translation to a Lagrangian frame. Alternatively, some frameworks integrate discrete agents that interact with continuum fields of nutrients, chemical signals, or mechanical stress, providing a flexible multiscale framework [18, 39]. Finally, in other hybrid approaches, continuous PDEs describe the global evolution of cell dis tributions, while discrete stochastic models capture rare switching events or lineage tracking [55, 22]. Correlating the equations and parameters of these two approaches is essential for building consistent, multiscale models of cell populations, and to facilitate the integration with data and between them. Many works such as [98, 99, 100, 101, 102, 103] derive a complete set of relationships between the parameters of a class of ABMs, incorporating cellular proliferation and death, as well as random and directed migration, with a class of continuum models that describe the combined effects of reaction, convection, and diffusion (see Figure 3 for an example of the equivalence between ABMs and continuum models in glioblastoma cultures in microfluidic devices).

Regarding the modeling of phenotypic plasticity, struc tured continuous epigenetic models have recently been derived from ABMs by taking the continuum limit, includ ing the dependence of migration and proliferation on the cell phenotypic state [105, 82, 77].

4.2. Structural Continuous Versus Internal Variables to

Model Epigenetic Transitions Discrete models concentrate each differential behavior in one specific phenotype. This may lead to a large amount of different phenotypes, or alternatively, if this number is maintained small, we may be forced to homogenize pheno typic differences that can be of interest in particular sit uations. This approach inherently assumes discontinuous boundaries between cell states, neglecting the continuous spectrum of transcriptional configurations that modern single-cell data reveal [7, 106].

To capture this gradual and history-dependent land scape without artificial discretization, frameworks must represent phenotypic evolution as a continuous process cou pled with environmental factors. This is typically achieved Fig 3: Equivalence between discrete and continuum mod els. Case example of the simulation of the formation of a necrotic core in glioblastoma (GBM) inside a microfluidic devices. Images reproduced from [104, 103, 42].

Through two main paradigms: extra-dimensional structural models [81, 25] and the internal variables approach [13, 14]. The primary advantage of extra-dimensional models lies in their completeness. They provide a detailed population level description of evolving heterogeneity, capturing full phenotypic distributions and explaining emergent phenom ena such as bimodality or selection-driven drift. This fa cilitates direct comparison with single-cell datasets, which increasingly report distributions rather than population averages. However, these models face significant hurdles; tracking individual cell trajectories —essential for model ing memory or hysteresis — is difficult unless further state variables are introduced. Furthermore, parameter identi fication is challenging, as the functional forms of vξ, Dξ, and R(c,x,ξ,t) are rarely directly measurable. Finally, they are computationally expensive, as the discretization of the extended domain is constrained by the curse of dimensionality as phenotypic dimensions increase. In contrast, the internal variables (IV) framework de scribes the mean evolution of the epigenetic state through a set of dependent variables and their evolution equations. Relative to the extra-dimensional family, the IV method offers several practical advantages: it naturally accommo dates long-term memory effects, aligns with typical experi mental regimes where low-order statistics are more readily available than full high-dimensional distributions, and remains modular. Adding new biological dimensions, such as a DNA-damage response state, augments the constitutive set without an exponential explosion in degrees of freedom, facilitating 3D tissue-scale studies coupled with nutrient diffusion and mechanics. The principal limitation of this method is its inability to resolve distributional shapes, such as tails or multimodality; when such features are critical, the model must be enriched with higher-order moments or reverted to a continuous distribution. Additionally, re ducing stochastic inheritance and division noise to scalar variables is a pragmatic but significant simplification of epigenetic variability.

The relationship between both approaches is clarified by interpreting the IV framework as a discretized projection of the structural model. If a continuous epigenetic coordinate follows a predefined probability density function (PDF), the set of moments defining that PDF can be identified as a finite set of internal variables whose evolution can be derived. This moment closure provides a natural bridge be tween both formalisms. Indeed, it is possible to define any finite-dimensional vector by means of a functional acting on a structured function f, λ = F(f) (note that statistical moments are particular cases). Under this setting, the evolution equation for λ is obtained immediately via the Gateaux derivative ˙ λ = D(f(·,t); ∂tf(·,t)). For instance, for models with the epigenetic state ξ as an independent variable and linear functionals (such as moments or collo cations), F[c](x,t) = Ωξ K(x,ξ)c(x,ξ,t)dξ, the Gateaux derivative is the functional itself. Thus, the evolution of the dependent variable ξ′ is given by In particular, for specific PDEs that have the closure prop erty for parametric families (e.g. exponential families), the evolution of f is strictly determined by the evolution of its defining parameters. This usually happens for linear evolution equations.

The computational efficiency of the IV approach is particularly evident in large-scale simulations: for a 3D hexahedral grid of N3 elements, a structural model with nξ dimensions (each discretized with m points) requires solving a system of N3 ×mnξ equations, whereas the IV approach solves only N3 spatial equations, followed by the efficient, parallelizable computation of (8 × nξ)N3 internal variables. Moreover, the Lagrangian formulation of the internal variable approach is usually simpler since knowing the initial distribution of epigenetic values at each spatial point is more complex. Conversely, integrating cell con centrations over the epigenetic coordinates and invoking suitable closures yields the evolution of the internal vari ables. This equivalence clarifies why internal variables often suffice in practice—especially when distributions remain narrow—and why extra-dimensional models become valu able when the distribution itself is the object of inference or control.

Figure 4 summarizes the main families discussed in this work to model epigenetic state, providing intuition to how they can be related to one another.

4.3. Eulerian Versus Lagrangian Formalisms to Model

Continuous Cell Processes Although the two continuous descriptions are equivalent mathematically, they emphasize different aspects. The anal ogy with continuum mechanics provides conceptual clarity. In fluid dynamics, the Eulerian formulation observes the flow field at fixed spatial locations, while the Lagrangian formulation follows individual fluid parcels. Thus, the Eu lerian approach captures the collective evolution of fields, enabling the analysis of global patterns, gradients, and fluxes. The Lagrangian formulation emphasizes the individ ual perspective, making it natural to model memory effects, and track the evolution of history-dependent internal vari ables. It is worth mentioning that in Eulerian approaches even if we take phenotypic plasticity into account, it is impossible to assign to a cell at a given RVE its history of perceived stimuli without solving before the transport equa tions, so a mechanistic approach for epigenetic transitions is harder. The following list of the differences between these two approximations highlights their main pros and cons. First, in the Lagrangian formulation, it is easier: 1. To ensure the conservation of quantities that affect a cellular RVE (e.g., concentration if there is no proliferation or death, and phenotype if there are no changes in it in a certain volume, although it may have changes per cell).

2. To track the evolution of a group of cells defined by any initial characteristic, such as their phenotype (tracking trajectories).

3. To track interfaces (e.g., cell groups with different initial phenotype).

4. To consider "mixing" problems (phenotypes close to each other but with different behavior).

5. To establish the relationship with discrete cell or dis crete phenotype models.

6. To consider memory, hysteresis or hereditary effects.

Fig 4: Summary of the different phenotypic plasticity models and how to derive equivalences between them. Red arrows indicate that the transition cannot be accomplished due to lack of information.

7. To know the variations in phenotype over time in relation to the phenotype history.

while in the Eulerian case it is easier:

1. To monitor the situation at a specific point or pheno type (e.g., how many cells of a given phenotype are present at a given moment, what is the flow of cells of a given phenotype at a specific spatial point).

2. To evaluate flow changes in a steady-state regime. 3. To set "wall" boundaries for movement and pheno type. 4. To solve problems where there are significant changes in position or phenotype, since the initial mesh cor responding to the phenotype domain may become highly distorted.

5. To study abrupt changes in the spatial or temporal evolution of the phenotype (phenotype waves).

The connection between these two formalisms is easy when knowing the movement x(X,t) and ξ = ξ(Ξ,t) and their inverses X = X(x,t) and Ξ = Ξ(ξ,t), which is not always possible to obtain in Eulerian formulations.

5. Conclusions

This review integrates biological and mathematical per spectives to elucidate how epigenetic transitions can be modeled across scales. We have traced how successive mod eling paradigms,from agent-based to continuum, seek to capture essential features of memory and adaptation. The analogy between physical and biological continua provides not only mathematical tools but also conceptual unity: the same notions of state variables, internal degrees of freedom, and transport apply to materials and to cells. Among the existing approaches for modeling pheno type as a continuum, which has been proven more bi ologically consistent, there are two main families. The extra-dimension formulation is powerful for exploring dis tributional structure, selection, and mixing. It embodies the field perspective, suitable for theoretical analysis of heterogeneity and for deriving macroscopic laws from mi croscopic mechanisms. The internal-variable formulation captures the same physics in a compressed, interpretable form suitable for practical simulation and integration with other models. Recognizing their equivalence and differ ences enables modelers to choose the appropriate level of description consciously rather than by convention.

Bridging the two — through theoretical analysis, hybrid modeling, and experimental calibration — will open the way toward a comprehensive, multiscale theory of epige netic dynamics capable of linking molecular mechanisms to population-level adaptation, and permit to better un derstand the actual differences between the models found in scientific literature.

Besides, our comparison highlights that the choice be tween Eulerian and Lagrangian frameworks is a matter of purpose rather than correctness. While the Eulerian description remains the standard for field evolution, a La grangian formulation becomes indispensable when tracking individual cellular histories and long-term memory effects. Future research should also focus on developing hybrid models that combine the strengths of both approaches. Extending internal-variable models to vectorial or tensorial internal spaces, incorporating stochasticity, and coupling them to gene-regulatory dynamics are important avenues for future work. Advances in single-cell multiomics, which can provide snapshots of the epigenetic distribution, will make it possible to estimate model parameters and val idate predictions quantitatively. Future extensions may include coupling to stochastic hybrid models, and with dimensionality-reduction techniques to infer the effective epigenetic coordinate from data, and data-driven identifi cation of drift and diffusion terms using physics-informed machine learning [107]. These developments will further strengthen the connection between mathematical modeling and experimental epigenetics, bringing predictive under standing of cellular adaptation closer to clinical application.

Acknowledgments

The authors would like to acknowledge the financial sup port from Beonchip S.L. through the project DIAMOOC IDI-20220819, the Spanish Research Agency (AEI), through the project PID2021-126051OB-C41 (GBM IMMUNE), the FORT23/00028 project of the ISCIII FORTALECE program, and the Government of Aragon, Spain (group DGA-T62_23R).

Declaration of Interest Statement The authors declare that they have no known competing financial interests or personal relationships that could have influenced the work reported in this paper.

Credit Autorship Contribution Statement Manuel Doblare Writing original draft, Writing review and editing, Conceptualization, Methodology, Inves tigation, Formal Analysis, Project administration, Funding acquisition. Marina Perez-Aliacar Writing- review and editing, Methodology, Investigation, Formal Analysis Jacobo Ayensa-Jiménez Writing- review and editing, Methodology, Investigation, Formal Analysis, Funding ac quisition.

Declararation of Ai-Based Tool Use During the preparation of this work Manuel Doblare used Google/-Gemini in order to help to translate an initial Spanish version to English. After using this tool, Manuel Doblare reviewed and edited the content as needed and takes full responsibility for the content of the publication. Besides, Jacobo Ayensa-Jiménez used Claude Sonnet 4.6 for reviewing the whole paper and the mathematical ex pressions and to homogenize the notation.