A cellular automaton model of osteogenic differentiation reveals identifiability limits of endpoint assays

1 Introduction

Biomineralization describes the formation and organization of inorganic phases in biological systems (Mann, 2001). In mammalian bone, this process is driven by osteogenic differentiation, a regulated program that has been studied extensively for decades (Stein et al., 1996; Izumiya et al., 2021). MC3T3-E1 mouse calvarial pre-osteoblasts are widely used as an in vitro model of osteogenesis because they mineralize consistently under defined culture conditions (Wang et al., 1999). Because MC3T3-E1 cells are lineage committed, variation in initial cell-state composition can alter differentiation trajectories and contribute to experimental variability (Croes et al., 2015).

Experimental characterization of osteogenic differentiation typically relies on endpoint assays that probe early, intermediate, and late stages of the program, including alkaline phosphatase (ALP) activity, collagen type I deposition, and mineralization measured by Alizarin Red S staining (ARS) (Stein et al., 1996; Franceschi et al., 1994; Ber et al., 2023; Rutkovskiy et al., 2016; Luttrell et al., 2019). These measurements are population-averaged and collected at discrete time points, and many assays are destructive, preventing repeated observation of the same cells over time (Tritschler et al., 2019). As a result, experimental data provide sparse and indirect observations of an underlying continuous differentiation process, limiting access to its dynamics and multi-scale interactions (Cook et al., 2024). In parallel, lineage commitment, matrix deposition, and mineral nucleation proceed asynchronously within MC3T3-E1 cultures (Fratzl-Zelman et al., 1998; Kou et al., 2026; Papantoniou et al., 2014).

Understanding differentiation dynamics requires linking cell-level processes such as proliferation, migration, and state transitions to experimentally observable outcomes. Cell-based computational models provide a framework for connecting individual behaviors to population-level readouts (Lee et al., 1995). Cellular automaton (CA) models represent multicellular systems using rule-based descriptions of cell processes and have been applied to growth and mineralization-related phenomena (Ghaffarizadeh et al., 2018; Truesdell and Saunders, 2020; Brüningk et al., 2019; Van Scoy et al., 2017). Existing in silico osteogenesis models typically reproduce temporal marker trajectories or spatial mineral patterns but do not explicitly examine whether standard assay panels contain sufficient information to distinguish alternative mechanistic hypotheses of lineage progression (Poorhemati and Komarova, 2024). Because these readouts provide only partial observations of the underlying process, the relationship between observable markers and differentiation dynamics constitutes a partially observed inverse problem (Jan et al., 2014). In such systems, multiple parameter combinations can yield similar observable behavior, leading to a many-to-one mapping between underlying dynamics and measured outputs (Gutenkunst et al., 2007).

Together, sparse endpoint sampling and many-to-one model behavior raise a central question: to what extent can standard endpoint measurements resolve underlying differentiation mechanisms? We frame this as a problem of partial observability in a stochastic cell-state system.

We study osteogenic differentiation under controlled media conditions. Standard osteogenic medium contains ascorbic acid (AA), which supports collagen synthesis, and β-glycerophosphate (β-GP), which provides phosphate for mineralization (Hwang and Horton, 2019). We introduce controlled perturbations using dexamethasone (Dex) and menaquinone-4 (MK-4), which act at partially distinct stages of differentiation: Dex promotes early lineage commitment (Jaiswal et al., 1997; Langenbach and Handschel, 2013), whereas MK-4 enhances osteogenic gene expression and matrix mineralization at later stages via γ-carboxylation–dependent pathways (Wu et al., 2019; Li et al., 2019). These stage-specific effects provide a basis to test whether standard assay panels can distinguish mechanistically distinct responses (Zhang et al., 2017).

We develop OsteoMin, a cellular automaton model linking latent MC3T3-E1 differentiation dynamics to assay readouts of ALP activity, collagen deposition, and mineralization. The model connects stochastic cell-state transitions, matrix accumulation, and mineral formation to experimentally accessible observables, enabling direct comparison with experimental measurements.

Using this framework, we examine how cell-state progression, spatial organization, matrix maturation, and media perturbations shape assay outcomes. Literature-parameterized simulations are compared with experimental endpoint data and complemented by module ablation and sensitivity analysis to identify mechanisms governing early and late readouts.

Finally, we assess practical identifiability under endpoint sampling by comparing endpoint features with time-resolved trajectories, testing whether standard assay panels can distinguish alternative perturbation families and quantifying how strongly measurements constrain underlying dynamics.

2 OsteoMin cellular automaton model

2.1 Overview of the OsteoMin model

OsteoMin was designed to establish a structural and interpretable link between latent osteogenic cell-state progression and commonly used endpoint assays. Rather than reproducing differentiation only at the level of generic growth dynamics or mineral pattern formation, the model directly connects stochastic transitions between pre-osteoblast and osteoblast states to experimentally measurable readouts. By aligning simulated outputs with assay observables, the framework enables mechanistic interpretation of endpoint measurements while also allowing assessment of how strongly such measurements constrain the underlying differentiation dynamics.

OsteoMin represents MC3T3-E1 monolayer cultures as a two-dimensional cellular automaton in which lattice sites are either empty or occupied by a single agent corresponding to a pre-osteoblast (pre-OB) or osteoblast (OB). Cells follow rule-based behaviors of migration, proliferation, differentiation, and apoptosis, while producing local extracellular fields corresponding to ALP activity, collagen matrix, and mineral (hydroxyapatite; HAp). These marker fields accumulate over time and are aggregated across the lattice to produce simulated assay readouts comparable to experimental ALP activity measurements, collagen immunofluorescence, and Alizarin Red S staining. The update scheme is illustrated in Figure 1.

A spatially explicit formulation is motivated by the observation that mineral deposition depends on local collagen density and osteoblast proximity (Fratzl-Zelman et al., 1998; Addison et al., 2015), processes that are not captured by population-averaged formulations. This spatial dependence is evident in the heterogeneous, nodular mineralization patterns observed experimentally (Figures 2h–k), consistent with locally clustered osteoblast activity (Mechiche A et al., 2016).

2.2 Spatial and temporal representation

The simulation domain is a square lattice of size n × n, with a default grid of 33 × 33 sites corresponding to an area of approximately 1 × 1 mm. Each lattice site represents the footprint of a single cell with spacing Δ x ≈ 30 μ m, consistent with the diameter of MC3T3-E1 cells. Periodic boundary conditions are applied: lattice sites at each edge are connected to those on the opposite edge, so that the domain has no preferred boundary direction or edge effects. Time is discretized as t k = k Δ t with timestep Δ t = 1 h, chosen to match reported migration timescales.

State updates follow a sequential asynchronous scheme. Cellular updates are applied in fixed order (movement, proliferation, commitment, apoptosis, marker production), while cells within each class are visited individually in randomized order at each timestep to minimize ordering bias.

Updates take effect immediately within each stage. During proliferation, target sites are reserved immediately, while daughter cells become active in the subsequent timestep. Cells undergoing apoptosis are removed before marker-field updates, ensuring extracellular quantities reflect the updated cell population.

Stochastic events were generated using NumPy’s pseudo-random number generator. Simulations used a fixed initial seed (12,345) unless otherwise stated. For experimental–simulation comparisons (Section 4.5) and latent initialization analyses (Section 4.3), independent randomly initialized seeds were used (n = 10 per condition).

2.3 Cell state dynamics

Pre-OBs are seeded at t 0 = 0 and can move, proliferate, and commit to the OB lineage. OBs remain stationary and deposit increased levels of assay markers. Apoptotic cells are removed from the lattice, leaving empty sites that can be occupied by pre-OBs.

2.3.1 Initialization

Simulations start by placing N 0 pre-OB cells on the lattice using random placement.

2.3.2 Migration

Pre-OBs perform a stochastic random walk on a square lattice ( Δ t = 1 h ), attempting a move to a randomly selected empty Moore-neighbor site with probability p m. Axial and diagonal directions are weighted to reduce lattice anisotropy and improve the rotational symmetry of long-time diffusion.

The migration probability p m was calibrated using migration data for MC3T3-E1 cells on flat substrates reported by Refaaq et al. (2020). The mean-squared displacement increases linearly with time, yielding an effective diffusion coefficient of ∼ 300–320 µm² h⁻¹. For a lattice spacing Δ x = 30µm, matching this hourly rate led to p m ≈ 0.27. Derivation of the migration mapping is provided in S1 Text, Supplementary Section 1.2.

2.3.3 Proliferation

Each pre-OB is assigned an individual division period τ i ∈ [ 38,52 ] h at birth (Cellosaurus, 2025). To represent cell-cycle progression, each cell carries an internal age variable a i k measuring the time since its last division. Cell age is updated at each timestep as a i k + 1 = a i k + Δ t.

When a i k ≥ τ i, the cell attempts to divide by placing a daughter pre-OB in a randomly selected empty Moore-neighbor site. Upon successful division, the parent cell age is reset to zero and the daughter cell is assigned a new division period τ i drawn independently. If no neighboring site is available, division is postponed and retried in subsequent timesteps.

2.3.4 Osteogenic commitment

Commitment from pre-OB to OB was modeled as a delayed stochastic transition. Following osteogenic induction at time t osteo, cells remain uncommitted for a maturation delay t d. After this delay ( t k ≥ t osteo + t d ), each eligible pre-OB transitions irreversibly to the OB state with probability p commit = 1 − exp ( − k commit Δ t ) per timestep, where k commit is the commitment rate.

The parameters t d and k commit were chosen to reproduce reported MC3T3-E1 differentiation kinetics, capturing the delay between osteogenic induction, early ALP activity, and subsequent mineralization (Dain Yazid et al., 2010; Izumiya et al., 2021). Additional implementation details are provided in S1 Text, Supplementary Section 1.3.

2.3.5 Osteoblast activity

OBs do not proliferate or move. They produce elevated levels of ALP, collagen, and mineral, which accumulate locally at the grid site.

2.3.6 Apoptosis

Apoptosis was modeled as a stochastic death process affecting pre-OBs. At each timestep, pre-OBs undergo apoptosis with probability p a = 1 × 1 0 − 3, representing low background turnover in culture (Bellows et al., 1998). Cells undergoing apoptosis are removed within the same timestep, while previously deposited matrix remains unchanged.

OB apoptosis was not included due to its low reported prevalence (0.6%) and limited quantitative constraints in vitro, which would introduce poorly identifiable parameters (Bellows et al., 1998).

2.4 Field dynamics and compound modulation

At each lattice site j, three state-associated quantities are tracked: ALP activity A j k, collagen matrix C j k, and mineral content H j k (representing HAp accumulation) at discrete times t k = k Δ t. These fields represent cumulative marker production associated with osteogenic differentiation and are aggregated across the lattice to obtain simulated assay readouts. All quantities are expressed in arbitrary units (AU), reflecting relative assay signal magnitudes.

ALP activity is primarily associated with cell-bound matrix vesicles. Marker dynamics are therefore treated as local, and diffusion and degradation are not explicitly modeled (Golub, 2011). Collagen accumulates locally as an extracellular matrix scaffold that supports mineral nucleation near OB regions (Franceschi et al., 1994; Fratzl-Zelman et al., 1998).

Ascorbic acid (AA) is modeled as a temporally varying medium component due to oxidative degradation (half-life ≈ 8 h) (Chepda et al., 2001), with degradation formulation and effective accumulation provided in S1 Text, Supplementary Section 1.4. β-GP was treated as constant due to its excess concentration (10 mM), regular medium exchange, and rapid diffusion relative to mineralization timescales.

Dexamethasone (Dex) and menaquinone-4 (MK-4) are treated as constant extracellular modulators refreshed during medium exchange. Their effects are represented using phenomenological compound–response functions chosen to reproduce literature-reported dose–response relationships within non-cytotoxic ranges (Dex 0–100nM; MK-4 0–10 µM), including intermediate optima observed for ALP activity, matrix production, and mineralization. Functional forms capture biphasic behavior while avoiding additional parameters that cannot be uniquely constrained from endpoint measurements (Langenbach and Handschel, 2013; Li et al., 2019; Alm et al., 2012; Kim et al., 2013; Della Bella et al., 2021).

2.4.1 Lattice state

Each lattice site j has staterepresenting empty sites, pre-OBs, and OBs, respectively. Binary occupancy variables indicate live-cell presence and OB identity,

S j k ∈ E m p , p r e − O B , O B ,

M j l i v e , k = 1 S j k ≠ E m p , M j O B , k = 1 S j k = O B .

Marker production rates are defined for each marker i ∈ { A , C , H } corresponding to ALP, collagen, and mineral fields. A basal production rate r b , i applies to all live cells, while OBs contribute an additional osteogenic production rate r o , i to the ALP, collagen, and mineral fields.

Osteogenic supplementation is introduced at time t osteo through the induction function defined in Equation 1:

u osteo t k = 1 t k ≥ t osteo (1)

which activates osteogenic production terms.

2.4.2 Compound modulation

Dex and MK-4 modulate osteogenic production through dimensionless response functions m i ( c D , c K ) acting on the osteogenic production rate r o , i, which is co-induced with the osteogenic medium. The compound modulation functions are defined below.

Both dexamethasone (Dex) and menaquinone-4 (MK-4) exhibit biphasic dose–response relationships in osteogenic cultures, with intermediate concentrations producing maximal effects on ALP activity, collagen deposition, and mineralization (Langenbach and Handschel, 2013; Wu et al., 2019; Li et al., 2019; Alm et al., 2012). Monotonic functions (Hill-type, logistic) cannot reproduce this behaviour without an additional inflection parameter. The log-parabolic formulation used here captures biphasic responses using a single shape parameter (b D , i or b K , i) per marker per compound, keeping the total number of constrained parameters minimal relative to the available calibration data. This is a phenomenological choice; other symmetric log-dose functions with similar shape properties would yield qualitatively equivalent results. The additive combination of Dex and MK-4 effects is the simplest assumption in the absence of data constraining interaction terms.

For each marker i ∈ { A , C , H } (ALP, collagen, HAp), Dex and MK-4 effects were modeled using biphasic log-dose response functions combined additively.

The marker-specific log-dose offsets are defined in Equation 2:where c D and c K denote Dex and MK-4 concentrations, and c D , 0 , i, c K , 0 , i are marker-specific reference concentrations associated with maximal response. For zero concentration, effectiveness was defined as g D , i = 0 and g K , i = 0.

x D , i = log 10 c D c D , 0 , i , x K , i = log 10 c K c K , 0 , i , (2)

The single-compound effectiveness functions are defined as clipped log-parabolic responses in Equation 3:such that g D , i , g K , i ∈ [ 0,1 ] represent relative compound effectiveness in log-dose space. Parameters b D , i , b K , i > 0 control the width of the biphasic response around the optimal concentration.

g D , i c D = max 0 , 1 − b D , i x D , i 2 , g K , i c K = max 0 , 1 − b K , i x K , i 2 , (3)

Combined modulation was assumed additive relative to the osteogenic baseline, as defined in Equation 4:where a D , i , a K , i > 0 specify the maximal contribution of Dex and MK-4 to osteogenic production for marker i. No explicit Dex–MK-4 interaction term was included.

m i c D , c K = max 0 , 1 + a D , i g D , i c D + a K , i g K , i c K , (4)

2.4.3 ALP dynamics

ALP activity evolves according to Equation 5:

A j k + 1 = A j k + Δ t M j l i v e , k r b , A + M j O B , k u osteo t k r o , A m A c D , c K . (5)

where r b , A denotes the basal ALP production rate in live cells, and r o , A denotes the additional osteogenic ALP production rate in OBs following osteogenic induction. The function m A ( c D , c K ) represents the dimensionless modulation factor describing the combined influence of Dex and MK-4 on ALP production.

The ALP field was initialized with baseline A j 0 = A 0, representing background media activity.

2.4.4 Collagen deposition

The osteogenic component of collagen production was modeled as AA-dependent using a Michaelis–Menten-type factor, as given in Equation 6:

C j k + 1 = C j k + Δ t M j l i v e , k r b , C + M j O B , k u osteo t k r o , C m C c D , c K c A A t k K A , C + c A A t k . (6)

The initial condition is C j 0 = 0. Here r b , C denotes the basal collagen production rate in live cells, and r o , C denotes the additional osteogenic collagen production rate in OBs after induction.

c A A ( t k ) denotes the AA concentration resulting from repeated medium supplementation and exponential decay, while K A , C defines the half-saturation constant for AA-dependent collagen production.

The modulation factor m C ( c D , c K ) describes the combined effects of Dex and MK-4 on osteogenic collagen production.

2.4.5 Mineralization dynamics

No basal mineral production is included ( r b , H = 0 ); mineral accumulation is exclusively osteogenic.

Mineralization occurs in OBs and depends on both extracellular matrix maturation and elapsed time since commitment. Commitment time is stored as a site-associated variable t j , c o m m i t for lattice sites occupied by OBs. The effective post-commitment time is defined in Equation 7:where t j , c o m m i t denotes the time at which the OB currently occupying site j committed. Delayed matrix maturation is represented by a mature collagen field C j , m a t k, which relaxes toward the deposited collagen level with timescale τ colaccording to Equation 8:

t eff , j k = t k − t j , c o m m i t . (7)

C j , m a t k + 1 = C j , m a t k + Δ t C j k − C j , m a t k τ col , C j , m a t 0 = 0 . (8)

Matrix-supported mineralization capacity follows the saturating dependence on matured collagen defined in Equation 9:Commitment-dependent mineral competence follows the Hill-type kinetics defined in Equation 10 as a function of elapsed post-commitment time:

g j , m a t k = C j , m a t k K col + C j , m a t k , (9)

h commit t eff , j k = t eff , j k n H τ H n H + t eff , j k n H . (10)

Mineralization evolves according to Equation 11:

H j k + 1 = H j k + Δ t M j O B , k u osteo t k r o , H m H c D , c K h commit t eff , j k g j , m a t k . (11)

Here H j 0 = H 0 denotes detectable background mineral signal. Parameter τ col controls collagen maturation, K col defines the half-saturation level for matrix-supported mineralization, and ( τ H , n H ) determine the timing and steepness of post-commitment mineral competence.

2.5 Literature-derived parameterization

Model structure and parameter values were constrained from prior mechanistic studies rather than inferred from the endpoint dataset presented here. Parameters were assembled from literature describing osteoblast proliferation, differentiation timing, matrix production kinetics, and Dex/MK-4 dose responses. Marker production rates and compound modulation terms were estimated by fitting reduced submodels to published ALP, collagen, mineralization, and dose–response data. Details of the fitting procedures and resulting parameter estimates are provided in S1 Text, Section 2.

Parameters used in simulations are summarized in Table 1. Probabilities are defined per timestep ( Δ t = 1 h ).

TABLE 1. Model parameters used in OsteoMin simulations.

A representative baseline simulation was performed over 28 days on a periodic 33 × 33 lattice with randomly seeded cells and no Dex or MK-4 perturbation. Given the assumed MC3T3-E1 cell footprint ( Δ x ≈ 30 μ m ), this domain corresponds to an approximate culture area of 1 × 1 mm, comparable to the spatial scale of a typical microscopy field of view. Cells undergo stochastic migration, pre-OB proliferation, apoptosis, and delayed osteogenic commitment following induction. ALP, collagen matrix, and mineral fields are produced locally and accumulate over time. Representative model outputs at assay-relevant times are shown for cell-state distribution (Figure 3a), ALP production (Figure 3b), and collagen–mineral co-localization (Figure 3c).

2.6 Statistical analysis

All experimental endpoint data are presented as individual points. Group differences were evaluated using the non-parametric Kruskal–Wallis test with two-sided Mann–Whitney U post hoc comparisons when significant. P-values were adjusted using the Benjamini–Hochberg procedure (false discovery rate control). Adjusted p-values are reported as q-values, with significance defined as q < 0.05 () and q < 0.01 (*).

Analyses were chosen according to the diagnostic objective of each model experiment. Concordance between simulated and experimental relative responses was evaluated using Kendall’s τ and directional agreement relative to the osteogenic condition. Separability of perturbation families was quantified using leave-one-out 1-nearest-neighbor classification of marker-specific trajectories (Euclidean distance), summarized by macro-F1. Parameter influence on assay-aligned endpoints was assessed using variance-based Sobol indices with Saltelli sampling (SALib).

2.7 Implementation and performance

The model was implemented in Python. Simulations ran on an AMD Ryzen 7 7435HS CPU (3.10 GHz), NVIDIA RTX 4060 Laptop GPU, and 16 GB RAM. Simulations used a 33 × 33 lattice (1,089 sites) seeded with 72 cells and run for 672 biological hours (28 days). Each 28-day simulation required approximately 22 s on a single CPU core (seed = 12,345).

3 Experimental methods

MC3T3-E1 pre-OB cells were cultured under basal and systematically varied osteogenic conditions to assess differentiation and mineralization responses.

3.1 Cell culture and osteogenic induction

MC3T3-E1 Subclone 4 cells were cultured in α-MEM (no ascorbic acid) supplemented with 10% FBS and 1% Antibiotic–Antimycotic at 37 °C and 5% CO 2. Cell concentration was determined by hemocytometer counting, yielding ∼ 4.1 × 1 0 5 cells/mL. For experiments, 12 mm glass coverslips (area ≈ 113 mm²) in 24-well plates were seeded at 8.0 × 1 0 3 cells/well (density ≈ 72 cells/m m 2) in 0.4 mL medium. Osteogenic supplements (ascorbic acid (AA), β -glycerophosphatee (β -GP), dexamethasone (Dex), menaquinone-4 (MK-4)) were added on Day 3, with medium changes every 48–72 h. The medium compositions are listed in Table 2. Detailed reagent information, including suppliers and catalog numbers, is provided in S1 Text, Section 6.

TABLE 2. Experimental media compositions.

All conditions contain ethanol carrier (0.10% v/v). Subscripts L and H denote low and high concentrations of dexamethasone (D; 10nM and 100 nM) and MK-4 (K; 0.1µM and 10 µM), respectively; subscript 0 indicates no supplementation.

3.2 Endpoint assays

Alkaline phosphatase (ALP) activity, collagen type I immunofluorescence, and Alizarin Red S (ARS) staining were performed on days 10, 14, and 28, respectively (n = 3 biological replicates). These assays and respective time points were selected as they represent complementary and sequential osteogenic differentiation stages in MC3T3-E1 cells. ALP activity is associated with early differentiation and typically peaks at days 7–10 in in vitro cultures. Collagen type I production was analysed at a time point when MC3T3-E1 cells are often reported to have established a mature extracellular matrix under osteogenic conditions, whie ARS staining was performed at a late time point when MC3T3-E1 cultures are known to exhibit robust calcium deposition and mineralized nodule formation (Yoon et al., 2025; Izumiya et al., 2021; Luttrell et al., 2019).

In addition to endpoint measurements, brightfield images were acquired using an Olympus IX81 inverted microscope with an sMAX04BM sCMOS camera and 4× and 10× objectives to monitor cell growth and mineralization throughout culture. Representative images are provided in S1 Text, Section 5.

3.2.1 Alkaline phosphatase activity

Intracellular ALP activity was assessed on day 10. Culture medium was aspirated, and cells were rinsed once with ice-cold Ca 2 +/M g 2 +-free DPBS to remove residual serum components. Cells were lysed using an ice-cold, phosphate-free buffer composed of 50 mM Tris, 1 mM MgCl 2, and 0.2% Triton X-100. Lysates were clarified by centrifugation at 10,000 × g for 10 min at 4 °C, aliquoted, and snap-frozen up to analysis time.

For enzymatic quantification, 25 µL of lysate was incubated with 75 µL of SIGMAFAST TM p-nitrophenyl phosphate (pNPP) substrate solution for 30 min at room temperature. The reaction was terminated by the addition of 25 µL of 3M NaOH, and absorbance was measured at 405 nm. Absorbance values were converted to ALP activity (mU⋅ mL − 1) assuming linear proportionality to enzymatic rate, using calf intestinal phosphatase as reference.

3.2.2 Collagen type I immunofluorescence

Extracellular matrix deposition and cellular organization were assessed by immunofluorescence staining of collagen type I, nuclei, and F-actin. Cells were rinsed with pre-warmed Ca 2 +/M g 2 +-free DPBS (pH 7.4) and fixed with freshly prepared 4% paraformaldehyde in DPBS for 15 min at room temperature, followed by three DPBS washes.

Cells were permeabilized with 0.1% Triton X-100 in DPBS for 5 min and blocked with 1% bovine serum albumin (BSA) in DPBS for 60 min at room temperature. Samples were incubated overnight at 4 °C with mouse anti-collagen type I primary antibody (1:300 in 1% BSA/DPBS).

After overnight incubation, samples were washed and incubated for 60 min at room temperature in the dark with goat anti-mouse Alexa Fluor 488 secondary antibody (1:500) and phalloidin–Alexa Fluor 568 (1:500). Nuclei were stained with DAPI (300 nM) for 15 min. Coverslips were mounted using Mowiol mounting medium and imaged the same day.

Images were acquired using a Leica SP8 confocal microscope equipped with an HC PL APO CS2 20×/0.75 IMM objective. Images of DAPI, actin, and collagen channels were acquired sequentially from 5 ROIs for all experimental conditions to allow quantitative comparison.

3.2.3 Alizarin Red S staining

Matrix mineralization was assessed on day 28 using ARS staining. Cultures in 24-well plates were rinsed with PBS and fixed with 4% paraformaldehyde in PBS for 15 min at room temperature. The ARS working solution was prepared at 1% (w/v), adjusted to pH 4.15 using 0.1M NaOH, and filtered through a 0.22 µm membrane to ensure selective binding to calcium-containing deposits. Fixed samples were incubated with the ARS solution for 15 min at room temperature under gentle agitation in the dark, followed by five washes with deionized water until complete removal of unbound dye, and then imaged.

Brightfield images were acquired using an Olympus IX-PHL inverted microscope equipped with a UPLFLN 4×/0.13 PH1 objective and an Olympus DP71 CCD camera. Nine fields of view were imaged per condition using identical exposure settings.

Mineralization was quantified as the fraction of ARS-positive pixels relative to total image area (Eggerschwiler et al., 2019). Images were processed using automated thresholding based on the complementary color ratio R / ( G + B ), where R, G, and B denote red, green, and blue channel intensities. Pixels with R / ( G + B ) ≥ 0.55 were classified as ARS-positive. Threshold sensitivity analysis is provided in S1 Text, Section 4.

4 Results

4.1 Experimental assay readouts

Osteogenic progression was evaluated at days 10, 14, and 28, corresponding to early enzymatic activity, extracellular matrix deposition, and late-stage mineralization.

4.1.1 ALP activity

ALP activity increased under osteogenic conditions relative to basal culture (Figure 2a). Dex and MK-4 supplementation produced variable responses without a consistent monotonic dose dependence. Differences between supplemented conditions were modest relative to replicate variability, and no statistically significant differences were observed.

4.1.2 Collagen deposition

Collagen type I area increased under osteogenic conditions compared to basal culture (Figure 2b). MK-4-containing conditions showed a modest increase relative to osteogenic medium alone, with substantial overlap between conditions.

Confocal imaging showed sparse collagen signal and less organized cytoskeletal structure in basal cultures (Figure 2d). Osteogenic conditions exhibited increased cell spreading and denser collagen networks (Figure 2e). Dex and MK-4 supplementation was associated with increased collagen signal intensity (Figures 2f,g) and organized actin filaments.

4.1.3 Mineralized area

Mineralized area fraction increased under osteogenic conditions relative to basal culture (Figure 2c). Differences between Dex/MK-4 conditions were moderate, with low Dex and MK-4 concentrations associated with increased mineralization.

ARS staining showed heterogeneous mineralization with localized nodular regions (Figures 2h–k). Spatial heterogeneity was further characterized using local variance analysis (S1 Text, Section 5) (Figure 3).

4.2 OsteoMin trajectories and module ablation

We first examined baseline temporal behavior and then quantified the contribution of individual mechanisms by selective module ablation. We use ablation in the computational sense: the model is rerun with a specific component disabled (e.g., k commit = 0, no apoptosis, no delay, or removal of the 2D spatial representation) and compared to the full model.

The reference spatial model generated assay-aligned trajectories for ALP, collagen, and mineralization (Figure 4a). ALP increased toward day 10, collagen increased toward day 14, and mineralization showed delayed onset with stronger growth toward day 28.

To assess mechanism-level contributions, each ablation was compared with the spatial reference model under identical baseline conditions. Nonspatial model details are provided in S1 Text, Section 3.2.1. Effects are reported as median fold change relative to the spatial reference model at assay endpoints across independent simulations, with significance from two-sided Mann–Whitney U tests and Benjamini–Hochberg correction (Figure 4b).

Module ablation showed distinct endpoint effects that clarify the functional role of individual mechanisms. Removing lineage commitment strongly reduced ALP (0.48×), collagen (0.59×), and mineral (0.27×), indicating that commitment is the primary source of osteoblast-dependent matrix production and mineral formation. Removing apoptosis increased all readouts (ALP 1.17×, collagen 1.11×, mineral 1.07×), indicating that turnover primarily limits cumulative matrix accumulation rather than regulating differentiation timing. Removing commitment delay produced minimal change in ALP (1.03×) and final cell number (1.00×), but modest reductions in collagen (0.92×) and mineral (0.97×), suggesting that differentiation timing influences matrix maturation windows more than early marker expression. Eliminating collagen-maturation delay left ALP and collagen unchanged (1.00×) while slightly increasing mineral (1.03×), consistent with mineral formation being constrained by matrix readiness rather than matrix abundance. Relative to the spatial reference, the nonspatial comparator increased early markers (ALP 1.80×, collagen 1.20×) and mineralization (1.17×), while modestly reducing final cell number (0.86×). This divergence indicates that spatial organization primarily affects late-stage mineral emergence: without spatial co-localization constraints, matrix accumulates rapidly but fails to achieve the local density and maturation coupling required for nucleation-like behavior. Together, these results indicate that commitment governs osteogenic capacity, apoptosis regulates cumulative production, and spatial exclusion with heterogeneous matrix maturation generates threshold-like mineralization dynamics that are not captured by population-averaged formulations.

4.3 Latent trajectories underlying similar endpoints

To examine how differences in initial cell organization affect model behavior, we simulated multiple initial seeding patterns on the same periodic (1 × 1 mm) lattice while keeping all biological parameters constant. All initialization patterns used the same total number of initial cells ( N 0 = 72 ), so differences reflect spatial organization rather than initial population size or global density. Across the initialization panel, seeding ranged from uniform random placement (C01) to structured motifs, including compact center clustering (C02), edge-biased and half-domain layouts (C03–C04), diagonal bands (C05–C06), corner clusters and multi-lobe patterns (C07–C10), lattice, stripe, and checkerboard geometries (C11–C13), radial spokes (C14), and a peripheral ring (C15) (Figure 5A).

For each configuration, ALP activity, collagen deposition, and mineralization were quantified at the corresponding assay endpoints and compared with a uniform random seeding control (C01). Mineral-field autocorrelation at 672 h was quantified using two-point correlation analysis (Figure 5B); details of the spatial autocorrelation metrics are provided in S1 Text, Section 4.

Relative to the reference condition (C01), endpoint variation was modest for ALP and mineralization but larger for collagen: ALP at 240 h ranged from + 1.1 % to − 9.5 %, collagen at 336 h from + 0.6 % to − 17.1 %, and mineral at 672 h from + 0.3 % to − 9.7 %. In particular, the central-cluster initialization (C02) showed the strongest reduction in mid-stage collagen accumulation. C10 slightly exceeded C01 because its quadrilobular seeding improved per-cell yield despite similar OB counts. Autocorrelation is shown as raw ξ: values span approximately 0.649 ± 0.010 to 1.382 ± 0.638, with C02 exhibiting the largest spatial correlation length. Overall, distinct spatial organizations can yield similar endpoint outcomes, indicating that endpoint assays underconstrain spatial aspects of differentiation and may obscure an important source of variability across experiments.

4.4 Effects of dex and MK-4 on osteogenic markers

Model responses to dexamethasone (Dex) and menaquinone-4 (MK-4) were evaluated across literature-reported concentration ranges (0–100 nM for Dex; 0–10 µM for MK-4) using single-factor parameter sweeps. Discrete dose grids are described in S1 Text, Section 3.2.

Endpoint values were extracted at assay-aligned time points (ALP: 240 h; collagen: 336 h; mineral: 672 h) and expressed as log 2 fold change relative to the osteogenic reference condition C2 (Dex = 0, MK-4 = 0).

To evaluate combined effects, Dex and MK-4 concentrations were varied jointly over a 10 × 10 grid spanning the same concentration ranges, yielding 100 conditions. For each condition, endpoint responses were computed using the same normalization procedure.

Dex single-factor responses exhibited biphasic behavior, with increases up to intermediate concentrations followed by mild declines at higher doses. ALP peaked at approximately 30nM, collagen at approximately 20nM, and mineral at approximately 60nM, with only marginal changes observed between 60 and 100 nM. MK-4 peak windows were ∼0.6–1 µM for ALP and mineralization (maximum at 0.6 µM), and 0.2–0.6 µM for collagen (maximum at 0.3 µM) (Figure 6a).

Combined responses were strongest within a bounded Dex–MK-4 region rather than at uniformly maximal single-factor doses (Figure 6b). In the two-dimensional mineral response surface (log 2 fold change vs. baseline osteogenic C2 condition), high responses form a broad plateau at intermediate-to-high Dex (∼20–100 nM) with low-to-moderate MK-4 (∼0.3–2 µM), while responses decline at high MK-4 (≳7–10 µM) across Dex levels. The local maximum is observed near Dex ∼40–60 nM and MK-4 ∼0.5–1 µM. Overall, the landscape indicates that multiple Dex–MK-4 pairs can achieve similarly high mineral outcomes within non-toxic concentration ranges.

4.5 Comparison of model and experimental endpoint responses

To evaluate agreement between independently obtained experimental measurements and literature-parameterized simulations, ALP, collagen, and mineralization endpoints were compared across the eight media conditions in Table 2. Model parameters were derived from prior studies and were not fitted to the experimental dataset presented here. The comparison therefore serves as an external consistency check of whether literature-supported differentiation dynamics produce qualitatively consistent endpoint trends under commonly used perturbations.

Responses were compared using a reference-centered effect metric relative to the osteogenic reference condition (C2) (Figure 7):

E = x − x C 2 | x | + | x C 2 | (12)

where x denotes the endpoint value under a given condition and x C 2 is the corresponding value under the reference osteogenic condition.

Concordance between simulated and experimental relative responses was evaluated using Kendall rank correlation ( τ ) and directional agreement (sign concordance) across non-reference conditions ( n = 7 ). Kendall τ quantifies consistency in condition ordering, while sign agreement measures concordance in response direction relative to the reference condition.

ALP showed limited directional agreement (14%) with weak rank concordance ( τ = 0.20 ). Collagen showed full directional agreement (100%) but negative rank concordance ( τ = − 0.14 ), indicating consistent response direction but poor ordering agreement across conditions. Mineralization (ARS) showed moderate directional agreement (43%) with moderate rank concordance ( τ = 0.43 ). Disagreement was concentrated in high-Dex conditions (C5–C7; Dex = 100 nM), especially for ALP and ARS, where simulated relative effects remained above C2 while experimental effects were reduced.

Overall, the comparison showed partial consistency but also substantial marker-specific discrepancies, especially under high-Dex conditions. This pattern supports the interpretation that endpoint agreement alone is insufficient to uniquely constrain the underlying mechanisms.

4.6 Endpoint observations weakly identify latent perturbation mechanisms

We evaluated whether staged endpoint assays contain sufficient information to distinguish mechanistic perturbation families arising from subclonal variation, spatial initialization, differentiation dynamics, and media-dependent modulation. These sources of variability are commonly associated with reproducibility differences in osteogenic experiments. Parameter ranges and sampling schemes for all perturbation families are provided in S1 Text, Supplementary Section 3.3.

We considered four perturbation families: spatial initialization, cell number, differentiation dynamics, and ascorbic-acid–mediated media effects, each comprising 50 parameterized conditions with 10 stochastic replicates. Replicate means were used to construct marker-specific endpoint and trajectory feature vectors. Values were standardized across conditions (for trajectories, per day), and Euclidean distances were used for leave-one-out nearest-neighbor classification of perturbation families. Performance was quantified using macro-F1, defined as the unweighted mean of per-family F1 scores.

Let z ( t ; θ , x 0 ) denote the latent differentiation trajectory determined by model parameters θ and spatial initialization x 0. Experimentally accessible measurements correspond to discrete observations of marker dynamics at assay-aligned timepoints (days 10, 14, and 28). For each marker m ∈ { A L P , C o l l a g e n , M i n e r a l } and condition i, we define (i) an endpoint feature and (ii) a truncated trajectory feature:with ( d ALP , d Collagen , d Mineral ) = ( 10,14,28 ).

x i , e n d m = m i d m x i , t r a j m = m i 0 , … , m i d m (13)

We emphasize that identifiability is assessed within the model class and observation scheme considered here. The results therefore quantify practical separability of mechanistic hypotheses under endpoint sampling, rather than intrinsic biological identifiability.

Macro-F1 ranges from 0 to 1, with a random baseline of 0.25 for this four-class problem. In the split-cell matrices (Figure 8), trajectory-based classifications show consistently stronger diagonal dominance than endpoint-based classifications, indicating improved recovery of true perturbation families. Endpoint scores remain modest (F 1 ALP@ 10 = 0.49, F 1 Collagen@ 14 = 0.40, F 1 Mineral@ 28 = 0.39), whereas trajectory-based scores are substantially higher (F 1 ALP ( 0 - 10 ) = 0.70, F 1 Collagen ( 0 - 14 ) = 0.86, F 1 Mineral ( 0 - 28 ) = 0.79).

Due to the sparsity of the available data and literature to estimate the parameters of Osteomin, these parameters have significant confidence intervals. Taking into account these intervals for the main parameters leads to a slightly reduced trajectory-based macro F 1 score (0.77), mainly driven by a reduction of the separability using ALP trajectories, while globally remaining higher than the endpoint macro-F 1 score (see S1. Supplementary Subsection 3.4).

This pattern reflects a collapse of distinct mechanistic perturbations onto similar endpoint observables, consistent with a many-to-one mapping from latent trajectories to measured outputs. Incorporating temporal information partially lifts this degeneracy, enabling improved discrimination between perturbation families.

4.7 Global sensitivity analysis identifies dominant regulatory parameters

Global sensitivity analysis was performed using variance-based Sobol indices implemented in SALib (Herman and Usher, 2017) to quantify parameter influence on simulated osteogenic endpoints. A total of D = 20 model parameter ranges were varied jointly across biologically plausible ranges (S1 Text, Section 3) using Saltelli’s extension of Sobol quasi-random sampling with base sample size N = 256, without computing second-order indices, resulting in N ( D + 2 ) = 5,632 model evaluations. The sampled parameter sets were drawn from uniform distributions over the specified bounds. Endpoints were extracted at assay-aligned timepoints (ALP 240 h, collagen 336 h, mineral 672 h).

Sensitivity results are summarized in Figure 9. Panels A–C show first-order and total-order indices for ALP, collagen, and mineral endpoints. Error bars indicate 95% confidence intervals.

Sensitivity analysis showed distinct parameter influence across the differentiation program. Early ALP levels were primarily influenced by proliferation timing, with the minimum division interval contributing the largest variance ( S T ≈ 0.62 ), consistent with ALP reflecting both lineage progression and population expansion (Luttrell et al., 2019; Izumiya et al., 2021). Collagen deposition depended on both division dynamics and basal and osteogenic collagen production, indicating combined effects of cell number and per-cell matrix secretion (Franceschi et al., 1994; Addison et al., 2015). Late mineralization was dominated by nonlinear mineral-specific parameters, particularly the Hill exponent governing mineral competence ( S T ≈ 0.41 ), consistent with delayed matrix maturation preceding mineralization (Addison et al., 2015; Zhu et al., 2024). The larger difference between first- and total-order indices for mineralization indicates interaction effects, suggesting that late-stage outputs depend on coupled upstream processes. Overall, different assay endpoints probe distinct components of the differentiation program.

4.8 Generalizability of the framework

Although OsteoMin was developed for osteogenic differentiation, the framework applies more broadly to in vitro systems in which cells undergo discrete state transitions, interact in space, and are primarily characterized through aggregate population-level measurements. In such settings, observable outputs reflect heterogeneous population behavior, while intermediate progression remains only partially resolved. Consequently, multiple latent dynamics may be consistent with similar endpoint readouts.

Rather than uniquely identifying mechanisms from endpoint data alone, the framework defines a constrained space of admissible latent dynamics. Given initial conditions and candidate transition rates, the model can be propagated forward to generate experimentally testable predictions across perturbations, doses, and sampling schemes. Comparison with independent measurements progressively excludes incompatible parameter regimes and supports hypothesis discrimination and experimental design in incompletely observed systems.

5 Discussion

OsteoMin treats ALP activity, collagen deposition, and mineralization as partial observations of an underlying stochastic cell-state process. Within this framework, endpoint assays correspond to an incomplete projection of differentiation dynamics, such that distinct mechanistic pathways—arising from different kinetic parameters or spatial initial conditions—can produce similar observable outcomes. This reflects a general limitation of inferring system dynamics from sparsely sampled, population-averaged measurements (Kitano, , 2002).

The observed non-identifiability has two contributing sources that should be distinguished. The first is structural: any model with partially compensatory processes—proliferation timing, commitment delay, and matrix maturation—will exhibit parameter non-identifiability regardless of the observation scheme, because multiple parameter combinations can produce similar population-level trajectories. The second is assay-dependent: endpoint sampling compresses trajectory information so that distinct latent dynamics can yield identical discrete readouts. The quantitative results presented here—macro-F1 ≈ 0.42 for endpoints and ≈ 0.78 for trajectories—are specific to the OsteoMin model structure, parameterization, and the perturbation families tested. However, the qualitative conclusion that temporal trajectories substantially outperform endpoint measurements in discriminating mechanistic hypotheses is expected to generalize broadly to any system with continuous heterogeneous dynamics assessed through sparse, population-averaged measurements.

This limitation is supported quantitatively by the model results. Using endpoint measurements alone, alternative mechanistic hypotheses were only weakly separable (macro-F1 ≈ 0.42), whereas incorporating temporal trajectories substantially improved discrimination (macro-F1 ≈ 0.78). This indicates that trajectory information, rather than terminal marker values alone, carries much of the mechanistic signal.

Spatial initialization further influences assay outcomes. The spatial formulation is motivated by the dependence of matrix accumulation and mineral nucleation on local cell density and neighborhood history, which are not directly represented in mean-field models without additional assumptions. Despite substantial variation in spatial organization, endpoint differences remained limited (<20% for collagen and ∼10% for mineral), indicating that spatial heterogeneity is only weakly reflected in standard assay readouts (Izumiya et al., 2021; Zhu et al., 2024). In the model, mineralization emerges from locally coupled matrix maturation and deposition, producing heterogeneous and clustered patterns consistent with experimental observations (Mechiche A et al., 2016). However, these differences are largely obscured when aggregated to endpoint measurements, such that distinct spatial organizations yield similar assay-level outcomes. The nonspatial formulation increased the mineral endpoint to 1.17× relative to the spatial reference, showing that spatial constraints influence the balance between cell growth and cumulative marker production. At the same time, the nonspatial formulation cannot capture local matrix‐mineral coupling or heterogeneous mineralization patterns. Establishing the biological importance of these spatial effects will require spatially resolved experiments, such as co-localization analysis of collagen and mineral at the single-nodule level, which were not included in the present experimental design.

Compound response analysis shows that the model reproduces biphasic dexamethasone responses, with intermediate concentrations associated with increased osteogenic markers (Alm et al., 2012; Chiang, 2025; Owen and Reilly, 2020). At the same time, multiple Dex/MK-4 combinations produce comparable endpoint responses, indicating that different regulatory regimes can map to similar observable states. This highlights that agreement with endpoint data alone does not uniquely constrain the underlying mechanism, even when model assumptions are consistent with known biology. The most notable model–experiment discrepancy occurs under high dexamethasone (100 nM, conditions C5–C7), where OsteoMin predicts ALP and mineralization above the osteogenic reference while experimental measurements show a reduction. The most likely biological explanations are dexamethasone-induced mitochondrial dysfunction and apoptosis at high concentrations (Kim et al., 2013; Ma et al., 2020), which reduce viable cell numbers and therefore ALP signal, and the sampling-time dependence of ALP activity: at high Dex, enzymatic activity may peak and decline before day 10, so that a single endpoint measurement underestimates the peak response (Owen et al., 1990). Neither mechanism is included in the current model: OB apoptosis was excluded due to limited quantitative in vitro constraints (Bellows et al., 1998), and the temporal profile of ALP under high Dex is not resolved by a single endpoint. Including these effects would require additional parameters that cannot be uniquely constrained from the available data, increasing model flexibility without resolving the underlying identifiability limitations (White et al., 2016). We therefore treat the high-Dex discrepancy as evidence of the model’s operating range: OsteoMin was designed to reproduce average-case osteogenic behaviour, and its quantitative predictions should not be extrapolated to near-cytotoxic concentration regimes.

Sensitivity analysis indicates that parameter influence shifts across the differentiation program. Early ALP variability is primarily associated with proliferation timing ( S T ≈ 0.62 ), whereas late mineralization depends on nonlinear parameters governing matrix competence ( S T ≈ 0.41 ). The larger difference between first- and total-order indices for mineralization suggests interaction effects, implying that late-stage outcomes depend on coupled upstream processes. As a result, endpoint measurements at a single time point provide limited ability to disentangle parameter contributions, further constraining practical identifiability.

Several limitations should be noted. OsteoMin was designed as a coarse-grained model to examine what can and cannot be inferred from endpoint assays, rather than as a fully biophysical description of osteogenic culture. The model is two-dimensional and does not include diffusion, matrix turnover, or mechanical feedback, all of which can influence osteogenesis in vitro (Golub, 2011; Engler et al., 2006). In this culture setting, however, diffusion is likely of secondary importance, as ascorbic acid is unstable in medium (half-life ∼8 h) (Chepda et al., 2001), and measured outputs primarily reflect cell-associated activity and locally retained matrix (Clarke Anderson, 2003). Dex and MK-4 were represented as phenomenological modulators of osteogenic activity, although both may also affect proliferation and lineage commitment through mechanisms not included here (Kim et al., 2013; Della Bella et al., 2021). Fitted parameters should therefore be interpreted as effective quantities rather than uniquely identifiable biological constants. A more detailed model would require substantially richer experimental data for calibration; without such constraints, added realism would mainly introduce unidentifiable parameters (White et al., 2016), consistent with the inability to capture inhibitory responses at high dexamethasone concentrations. The log-parabolic dose–response functions and the additive combination of compound effects are phenomenological choices that adequately reproduce the calibration data but are not derived from mechanistic first principles; alternative parameterizations within the same functional class would yield qualitatively similar predictions. The current comparison uses a single literature-derived parameter set and does not formally propagate parameter uncertainty into prediction intervals. The simulated trajectories therefore represent the model’s average-case behaviour rather than a distribution of plausible outcomes. OsteoMin is structurally compatible with ensemble-based uncertainty propagation—for example, by sampling from parameter distributions informed by experimental variability in division periods, commitment rates, and marker production rates—but doing so rigorously requires a more complete quantification of biological variability than is currently available from the literature. We regard this as an important direction for future work; richer experimental datasets with quantified within-replicate distributions would also improve the ability to distinguish alternative differentiation mechanisms.

Relative to earlier osteogenesis models, OsteoMin explicitly treats experimental assays as partial observations of a stochastic cell-state process. Previous cellular automaton models reproduced spatial mineralization patterns without systematically linking intermediate markers to assay-level observables (Van Scoy et al., 2017), while continuum and ODE approaches capture aspects of matrix maturation but typically assume homogeneous populations (Komarova et al., 2015).

These results indicate that most identifiable information lies in marker dynamics rather than endpoint values, which compress the underlying progression and limit mechanistic resolution. Viewed as a partially observed inverse problem, endpoint assays do not identify the underlying dynamics but instead define a set of admissible trajectories. The role of the model is therefore not to recover a unique mechanism, but to map this space of consistent dynamics and test how different hypotheses or perturbations project onto observable endpoints. Within this framework, separability depends on when and what is measured, indicating that experimental design—particularly assay timing—can be optimized to better resolve otherwise indistinguishable mechanisms.

The identifiability analysis points to concrete experimental design implications. The trajectory-based classification results (Section 4.6) indicate which markers and time windows carry the most discriminative information. ALP trajectories over days 0–10 ( F 1 ALP ( 0 - 10 ) = 0.70 ) substantially outperform the single day-10 endpoint ( F 1 ALP@ 10 = 0.49 ), suggesting that two or three intermediate ALP measurements (e.g., days 3 and 7) would improve discrimination with minimal additional cost. Collagen time-courses show the highest separability of any individual marker ( F 1 Col ( 0 - 14 ) = 0.86 ), making longitudinal collagen quantification particularly valuable. Additional markers that could further improve separation include total cell density (which separates proliferation-driven from differentiation-driven trajectories) and spatially resolved mineral quantification (e.g., autocorrelation length from ARS images). OsteoMin can be used directly to simulate expected classification performance under alternative measurement panels before experiments are performed, enabling in silico optimization of assay design.

6 Conclusion

We developed OsteoMin, an interpretable cellular automaton that structurally links stochastic osteogenic cell-state progression to commonly used endpoint assays of ALP activity, collagen deposition, and mineralization. By directly connecting latent differentiation dynamics to assay-level observables, the framework enables mechanistic interpretation of population-averaged readouts within a spatially resolved model of lineage progression.

Model results indicate that similar endpoint measurements can arise from distinct underlying trajectories, showing that standard assay panels constrain broad progression trends but often do not uniquely identify differentiation mechanisms. OsteoMin therefore provides both an interpretable model of osteogenic differentiation and a quantitative basis for assessing practical identifiability limits in endpoint-driven experiments, supporting the design of measurements that better resolve latent cell-state dynamics.

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/fsysb.2026.1872774/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Kevin Leder, University of Minnesota Twin Cities, United States

Reviewed by: Javiera Cortés-Ríos, Pontificia Universidad Católica de Chile, Chile

Carlo Galli, University of Parma, Italy