# Stiffness-sensitive gene regulation in human mesenchymal stem cells: Modelling mechanotransduction to predict mineralization and bone protein expression

Jean-Philippe Berteau, Abdennasser Chekroun, Laurent Pujo-Menjouet, Kevin Yueh-Hsun Yang

► To cite this version:

Jean-Philippe Berteau, Abdennasser Chekroun, Laurent Pujo-Menjouet, Kevin Yueh-Hsun Yang. Stiffness-sensitive gene regulation in human mesenchymal stem cells: Modelling mechanotransduction to predict mineralization and bone protein expression. Journal of Theoretical Biology, 2026, 618, pp.112284. 10.1016/j.jtbi.2025.112284. hal-05376493

HAL Id: hal-05376493

https://hal.science/hal-05376493v1

Submitted on 21 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Stiffness-sensitive gene regulation in human mesenchymal stem cells: modelling mechanotransduction to predict mineralization and bone protein expression

Jean-Philippe Berteau<sup>a,b,c</sup>, Abdennasser Chekroun<sup>d</sup>, Laurent Pujo-Menjouet<sup>e,∗</sup>
and Kevin Yueh-Hsun Yang<sup>f</sup>

November 21, 2025

<sup>a</sup> Department of Physical Therapy, College of Staten Island, City University of New York, 2800 Victory Blvd, Staten Island, 10314, NY, USA.

<sup>b</sup> New York Center for Biomedical Engineering, City College of New York, City University of New York, 160 Convent Ave, New York, 10031, NY, USA.

<sup>c</sup> Nanoscience Initiative, Advance Science Research Center, City University of New York, 85 St Nicholas Terrace, New York, 10031, NY, USA.

<sup>d</sup> Laboratoire d'Analyse Nonlinéaire et Mathématiques Appliquées, University of Tlemcen, Tlemcen 13000, Algeria.

<sup>e</sup> Université Claude Bernard Lyon 1, CNRS, Centrale Lyon, INSA Lyon, Université Jean Monnet, I CJ UMR5208, Inria, 69622 Villeurbanne, France.

<sup>f</sup> Translational Tissue Engineering Center, Johns Hopkins University School of Medicine, Baltimore, 21213, USA.

## Abstract

The goal of our study was to establish how a specific part of the bone Gene Regulatory Network (GRN) controls mineralization in response to stiffness. We hypothesized that a system of differential equations model stiffness-sensitive gene regulation in human mesenchymal stem cells through the epistatic genetic interactions between stiffness (e.g. WNT-β catenin pathway) and five of the main transcription factors and bone proteins (e.g. RUNX2, BSP, OSX, OC, and OPN). To test this hypothesis, we (i) performed in-vitro experiments culturing bone cells on different stiffness, (ii) adapted our previously published model from being continuously time-dependent to continuously stiffness-sensitive, and (iii) simulated protein production in function of stiffness and other protein production from the best estimate of parameters coming from the experimental work. Our experimental findings reveal a non-parametric relationship between stiffness and RUNX2 production, with no discernible linear trends for other proteins. Modeling results demonstrate that continuous variations in stiffness enable simulation of bone GRN gene expression, fitting our novel experimental dataset. Specifically, our computational results indicate that OPN production peaks at low stiffness (8 kPa), while RUNX2, OSX, and OC achieve maximum production at higher stiffness levels (64 kPa). This alignment underscores the model's capacity to replicate experimental data accurately. Additionally, our approach predicts that WNT-β-catenin activation serves as an enhancer for OPN and BSP production. The
model also highlights a negative feedback-like interaction between OC and BSP production. Stiffness variations were shown to have a significant impact on OC and BSP production and a moderate effect on OPN production. By employing a stiffness-sensitive gene regulation model, we provide insights into one of the mineralization patterns through the prediction of bone protein expression dynamics.

**Keywords:** bone modeling, Mathematical model, Ordinary differential equations, gene regulation, Mineralization.

## Introduction

To build bone structure - mineralizing a collagen matrix and remodelling it over time -, bone cells initiate a self-assembly process monitored by the bone genetic regulatory network (GRN). Evidence has proposed several theoretical approaches to understand what controls bone mineralization and remodelling [1, 2, 3, 4]; unfortunately, few of them address the main GRN components, their behaviour, and the order they act on each other [5, 6, 7]. In our previous work, we detailed the role played by each production of some transcription factors, enhancers, and inhibitors of mineralization initiated by mechanical forces acting on bone cells - called the canonical pathway (Wingless/Beta Catenin, e.g. WNT-β catenin pathway, mechanotransduction) - to build a bone [8]. However, the model proposed is limited to a theoretical approach, and some parts of the bone GRN are missing - such as one of the major enhancers of mineralization called BSP (bone Sialoprotein).

Regarding interactions between stiffness and gene expression, the canonical mechanotransduction pathway (e.g. WNT-β catenin pathway) regulates how bone cell mechanoreceptors stimulation initiates bone mineralization [10, 11, 12, 13, 14]. It can be activated at different scales by (i) the stiffness of the extracellular bone matrix [15], (ii) shear forces coming from fluid movements surrounding the cell [16], (iii) vibration propagated along the bone diaphysis [4, 17], (iv) tension coming from muscles' contractions [18], or (v) compression coming from body movements acting under gravity [19]. Here, we are interested in the first mechanism of action which is how alteration of the bone matrix stiffness activates the bone mineralization pathway. Regarding the bone mineralization pathway, a sequence of genes are known to be essential: (i) Runt-related transcription factor 2 (RUNX2), (ii) Osterix (OSX) (iii) Special AT-rich sequence-binding protein 2 (SATB2) which in turn initiates the production of iv) two enhancers of mineralization, namely, bone sialoprotein (BSP) and alkaline phosphatase (ALP), and of v) two inhibitors of bone mineralization called osteocalcin (OC) and osteoponin (OPN) [20, 21, 22, 10].

Mathematical modelling of networks emerges as a crucial complementary approach, capable of integrating qualitative and quantitative information about network architecture and parameters to offer an integrative platform for comprehending the outcomes of diverse genetic perturbations and generating novel predictions. Evidence has already proposed several theoretical approaches to bone mineralization and remodelling [1, 2, 3]), but few of them [8, 9] describe the role at the level of the cell. In our previous paper [8], we detailed each production of transcription factors, enhancers, and inhibitors of mineralization by using the canonical pathway of bone GRN activation (Wingless/Beta Catenin) and modelled transcription factors and bone protein production through Michaelis-Menten and Hill function. As a result, we proposed two nonlinear differential equations that fit well with the previously published experimental data. Compared to empirical evidence - coming from in-vitro experiments culturing bone cells on different stiffness and over time [26] -, the two best systems used factors of inhibition from the start of the activation of each gene. The difference between the two systems lay in the BSP equation and two ways to activate and reduce its production and revealed negative indirect interactions from negative feedback loops and the recently depicted micro-RNAs, providing the first theoretical evidence of osteoblast self-inhibition after activation of the genetic regulatory network controlling mineralization. However, the model lacked BSP data - a key contrib-
utor to mineralization [27, 28, 29] - and the time dependence of our dataset appeared to be one of the significant limitations. Indeed, evidence has shown that stiffness alterations of the extracellular matrix from increased collagen crosslinking and mineralization is the primary trigger of bone GRN activation [15, 30, 31] and not time. Thus, stiffness should be the primary input in a theoretical model that addresses epistatic genetic interactions that govern bone mineralization through mechanotransduction and the limits of our previous approach [8] was that stiffness was not included in the model, explicitly.

The goal of our study was to establish how a specific part of the bone Gene Regulatory Network (GRN) controls mineralization in response to stiffness input. We hypothesized that a system of differential equations model stiffness-sensitive gene regulation in human mesenchymal stem cells through the epistatic genetic interactions between stiffness (e.g. WNT-β catenin pathway) and five of the main transcription factors and bone proteins (e.g. RUNX2, BSP, OSX, OC, and OPN). Our hypothesis was that these interactions follow a system of nonlinear differential equations. To test this hypothesis, we (i) performed *in-vitro* experiments culturing bone cells on different stiffness, (i) adapted our previously published model [8] from being continuously time-dependent to continuously stiffness-sensitive, and (iii) simulated protein production in function of stiffness and other protein production from the best estimate of parameters coming from the experimental work.

## Materials and methods

We used the following three-step method:

* Biological Experimentation: we investigated transcription or protein production using a standard quantitative real-time PCR procedure after Human Mesenchymal stem cells (hMSCs) culture under four different stiffness: 8, 16, 32 and 64 kPa.

* Mathematical Modelling: we adapted our previously published model from being continuously time-dependent to continuously stiffness-sensitive.

* Simulation: we simulated bone GRN expression based on stiffness input (i.e., mechanotransduction) and investigated protein-protein interactions. In this study, we utilize the term “Gene Regulatory Network (GRN)” to refer to a model representing the intricate interactions among genes and their products that control gene expression levels within a specific biological system, focusing on bone. The GRN is a valuable tool for elucidating regulatory relationships between genes and proteins, including transcription factors, mineralization enhancers, and inhibitors, as they pertain to gene expression in the bone GRN. This graphical representation aids in our understanding of various biological processes, such as development, responses to stimuli, and diseases, with a primary emphasis on bone mineralization dynamics.

## Biological Experimentation

This experiment aimed to study the differentiation of Human Mesenchymal stem cells (hMSCs) into osteoblast. To do so, we used a qRT-PCR standard procedure to determine the fold change of proteins and transcription factors of hMSCs cultured on a petri dish coated with different stiffness and compare it to the fold change in hMSCs not cultured (control). qRT-PCR measures RNA transcripts rather than proteins or transcriptors factors themselves; therefore, in our study of the bone genetic regulatory network, we used RNA levels as proxies for protein expression [41]. Four coatings with different stiffness values (8, 16, 32, and 64 kPa, a sequence that increases by a factor of two between each value, covering the commonly used minimum and maximum stiffness thresholds for
<figure><img src="image_2.png" /><figcaption>Figure 1: Representation of the WNT-β catenin pathway used to trigger the epistatic genetic interactions modelled here that govern bone mineralization through mechanotransduction (transcription factors RUNX2, OSX, and bone proteins BSP, OPN, OC) with their respective numbered chromosomes - BSP and OPN are on the same (4). The red arrow signifies the downstream effect of Wnt activation, where nuclear β-catenin directly influences transcription within the bone GRN. Based on the research results in [32, 26], the figure illustrates how mineralization increases the stiffness and the associated feedback loop of the extracellular matrix on the bone Gene Regulatory Network (GRN). The gene regulation network and its interactions are described in grayscale with an example of trend behavior for each protein when stiffness increases (RUNX2 increased, OSX decreased, BSP decreased, OPN increased, OC decreased)</figcaption></figure>

mechanotransduction) were used to trigger transcription factors and protein production (RUNX2, BSP, OSX, OC, and OPN). Three replicates were performed, and we determined the relative fold changes in bone-specific proteins and transcription factors at a timepoint of *t* = 336 hours (= two weeks).

Sample preparation and Cell culture

Human mesenchymal stem cells (hMSCs) were purchased from ATCC (Manassas, VA, USA). The hMSCs at passage four and seeding density of 1428 cm² were cultured in α-MEM medium (Fisher) containing 16.7% FBS (Foetal Bovine Serum) (ThermoFisher, Waltham, MA, USA), one per cent PSF (penicillin/streptomycin) (ThermoFisher, Waltham, MA, USA) and one per cent L-glutamine (ThermoFisher, Waltham, MA, USA) in a 175 cm² cell culture flasks. The hMSCs culture was
maintained at 37°C in a humidified CO₂ incubator for 10 days with medium renewal every day for the first two days and every three days for the rest of the culture. The hMSCs were seeded with DMEM high glucose medium (ThermoFisher, Waltham, MA, USA) containing 10% FBS (ThermoFisher, Waltham, MA, USA), one per cent PSF (ThermoFisher, Waltham, MA, USA), one per cent L-glutamine (ThermoFisher, Waltham, MA, USA). 4×10⁵ cells were used as control samples. We used CytoSoft® 6-well plate contains a consistent 0.5 mm thick layer of silicone gel in each well, with its elastic modulus: 8, 16, 32 to 64 kPa to offer a range of substrate rigidities for cell culture, allowing for elasticity modulation (i.e. changing the stiffness of the all substrate). Since previous results [26] suggest that culturing hMSCs on substrates with a stiffness of 62-68 kPa promotes their differentiation into osteoblasts, here the chosen stiffness values (8 kPa, 16 kPa, 32 kPa, and 64 kPa) represent a range of mechanical properties that bone cells might encounter in various physiological conditions during growth, aging, mineralization. These values, while much lower than the GPa range typically reported for macroscopic bone tissue stiffness, are designed to mimic localized microenvironments or specific experimental substrates. This controlled setup allows us to test the hypothesis that stiffness-sensitive gene regulation influences mineralization.

### qRT-PCR sample

RNA purification and quantification were performed using TRIzol and GlycoBlue (ThermoFisher, Waltham, MA, USA). CDNA synthesis followed SuperScript™ III First-Strand Synthesis SuperMix for qRT-PCR protocol's directions from Invitrogen (Carlsbad, CA, USA). CDNA dilution and qRT-PCR followed the TaqMan Gene Expression assay protocol's directions from Applied Biosystems (Foster City, CA, USA).

### Acquisition protocol

The differentiation of the hMSCs into osteoblasts was assessed by measuring the relative fold change in bone-specific proteins and transcription factors, RUNX2, BSP, OSX, OC, and OPN, using a standard quantitative real-time PCR procedure. We used the following probes RUNX2 (Hs01047973_m1) OC/BGLAP (Hs01587814_g1) IBSP/BSP (Hs00913377_m1) SPP1/OPN (Hs00959010_m1) Osterix/SP7 (Hs01866874_s1) GAPDH (Hs02786624_g1) through TaqMan™ (ThermoFisher, Waltham, MA, USA).

### Data extraction

PCR results were extracted using the One Step software v from Applied Biosystems (Foster City, CA, USA) under a table format. The qPCR results were determined based on the threshold cycle numbers (CT) detected in the instrument and the comparative CT method was used to calculate the fold change of each gene investigated relative to the control.

### Statistical analysis

PCR data were analysed for (i) normality using Shapiro-Wilk test; (ii) parametric or non parametric relationship using Pearson and Spearman correlation indices, respectively. We used IBM Corp. Released 2023. IBM SPSS Statistics for Windows, Version 29.0.2.0 Armonk, NY: IBM Corp.

### Mathematical Modelling

We set up a new model - where each variable and parameter are described in Table 1 - to assess GRN protein concentrations with respect to several stiffness points but only a one-time step. We introduced a simplified model for the stiffness feedback loop (Fig. 1). Thus, we split our previous system [8] into five sub-models of two equations describing the increase of stiffness and the production
of each transcription factor, nucleator, and inhibitor. Thus, we built a novel stiffness-sensitive model with four stiffness values (8, 16, 32, and 64 kPa) for simulating the impact of stiffness on transcription or protein production (RUNX2, BSP, OSX, OC, and OPN). The objectives of the proposed model are then to (i) reject or accept the interaction assumptions between a continuous stiffness and the different GRN elements, (ii) fit the data according to stiffness by eliminating the role of time for RUNX2, BSP, OPN, OSX, and OC. Each model equation is structured in two parts [8]: a positive part representing the growth of stiffness and the activation of polypeptides, and a negative part describing the decrease of stiffness impact the concentration loss of polypeptides (by degradation or recruitment). Here, we neglect the stiffness action in the degradation part to homogenise the sub-equations.

<table><caption>Table 1: Non-negative variables and parameters considered in our models. Note that all variables $E_y$, $R_U$, $B_S$, $O_P$, $O_C$ and $O_X$ are time-varying for $t \ge 0$.</caption><thead><tr><th>Variables</th><th>Description</th></tr></thead><tbody><tr><td>$E_y$</td><td>Stiffness</td></tr><tr><td>$R_U$</td><td>RUNX2 concentration</td></tr><tr><td>$B_S$</td><td>BSP concentration</td></tr><tr><td>$O_P$</td><td>OPN concentration</td></tr><tr><td>$O_C$</td><td>OC concentration</td></tr><tr><td>$O_X$</td><td>OSX concentration</td></tr><tr><th>Parameters</th><th>Description</th></tr><tr><td>$n_S, n_{R_U}, n_{B_S}, n_{O_P}, n_{O_C}, n_{O_X}$</td><td>Sensitivity constants of Hill functions</td></tr><tr><td>$\mu_S, \mu_{R_U}, \mu_{B_S}, \mu_{O_P}, \mu_{O_C}, \mu_{O_X}$</td><td>Degradation rate or interaction through a non-canonical pathway with other components</td></tr><tr><td>$k_S, k_{R_U}, k_{B_S}, k_{O_P}, k_{O_C}, k_{O_X}$</td><td>Saturation levels for each concentration</td></tr><tr><td>$\alpha_R, \alpha_B, \alpha_C, \alpha_P, \alpha_X$</td><td>Saturating weight</td></tr></tbody></table>

## Positive part and general model

For the positive part, we describe the stiffness saturation and the polypeptide concentration by a version of the Hill function, inspired by the one in [8]. The expression of the function is, for all $x \ge 0$ and $y \ge 0$,

$$f_x(y) := \frac{ky^n}{b^n + x^n y^n}, \quad k > 0, b > 0, n \in (-\infty, +\infty), \qquad (1)$$

where $y$ is the stiffness or concentration of polypeptide and $x$ is the concentration of polypeptide interacting with $y$. The parameters $k$ and $b$ are real positive numbers, and $n$ is a real number describing the saturation level, activation coefficient, and sensitivity term (see [8] for the details). When $n$ is positive, the function describes the increase in the first phase, and when $y$ is equal to the activation coefficient $b$, the function reaches the half maximum point, and then the function saturates in the second phase. We add the action of polypeptide $x$ or the action of stiffness to the function of Hill by the product $x^n y^n$, which slows down or accelerates the function according to the reaction of $x$ at the first phase. When $n$ is negative, the reactions are the inverse of the previous case, that is, functions are decreasing with respect to $x$ or $y$ (see subsection for its explicit formulation. Fig. 2 ($n$ positive) and Fig. 3 ($n$ negative) show the influence of Hill coefficient $n$ and the interaction of the element $x$ in the Hill function. At this stage, it is crucial to emphasize that Hill functions are commonly used in biological models to represent saturation effects and molecular interactions [33]. Given this widespread application, it is logical to adopt a Hill function to describe polypeptide interactions in
our model. Specifically, we draw on Ferrell's research [33] to justify our choice, particularly in the context of chemical kinetic cascades that require reaching a certain concentration threshold before activation occurs.

While alternative S-shaped functions exist, we view them as primarily phenomenological and not biologically grounded. Functions such as the arctangent, hyperbolic tangent, or various polynomial or exponential ratios might have been considered, but they may lack the mechanistic justification that the Hill function provides from our point of view. Furthermore, given the limited understanding of the biological processes involved and the sparse experimental data available — both in terms of time and stiffness — our approach is to minimize speculative assumptions and instead present the simplest, most mechanistically plausible model. Thus, we deliberately choose the Hill function due to its strong alignment with biological mechanisms.

Regarding the influence of $E_y$ or $x_i$ on themselves, we propose a biological hypothesis where self-regulation is framed as a negative feedback loop. Specifically, we suggest that increasing stiffness affects the growth rate of the $x_i$ component, ultimately leading to a decelerating effect. This interaction is central to our conceptual framework.

Finally, the term $y^n$ in $x^n y^n$ is a standard feature of the Hill function. The novelty in our approach lies in replacing a constant with $x^n$, a decision informed by our previous work [8], which identifies a deceleration effect tied to different polypeptides.

<figure><img src="image_3.png" /><figcaption>Figure 2: The case where $n$ is positive: Hill function $y \mapsto f_x(y) = \frac{ky^n}{b^n+x^n y^n}$ with $b = 2$ and $k = 1$. In this case, the function $f_x(y)$ is increasing with respect to $y$. Left: influence of $n$ (positive) with $x = 1$; the saturation curve becomes steeper when Hill's coefficient increases. Right: influence of the element $x$ with $n = 1 > 0$; when $x$ increases, the function reaches saturation faster.</figcaption></figure>

The general system to describe the interactions between the stiffness and each element of GRN is then given by the following system, for $t \ge 0$,

$$ 
\left\{ 
\begin{aligned}
\frac{dE_y}{dt} &= \frac{k_S E_y^{n_S}}{b_S^{n_S} + \sum_{i=1}^{5} (\alpha_i x_i^{n_S}) E_y^{n_S}} E_y - \mu_S E_y, \\
\frac{dx_i}{dt} &= \frac{k_{x_i} x_i^{n_{x_i}}}{b_{x_i}^{n_{x_i}} + E_y^{n_{x_i}} x_i^{n_{x_i}}} x_i - \mu_{x_i} x_i, \quad i = 1, \dots, 5,
\end{aligned}
\right.
\quad (2) $$

where $E_y$ represents the stiffness and $x_i$ represents for simplicity the different elements of GRN (RUNX2, BSP, OPN, OSX and OC). Note that for clarity and simplicity we write $E_y$ and $x_i$ instead of $E_y(t)$ and $x_i(t)$. This simplification of the notation is applied throughout the rest of the paper.
For $t \ge 0$, we thus denote

$$x_1(t) := R_U(t), \quad x_2(t) := B_S(t), \quad x_3(t) := O_C(t), \quad x_4(t) := O_P(t) \quad \text{and} \quad x_5(t) := O_X(t).$$

Also, we set

$$\alpha_1 := \alpha_R, \quad \alpha_2 := \alpha_B, \quad \alpha_3 := \alpha_C,$$

$$\alpha_4 := \alpha_P \quad \text{and} \quad \alpha_5 := \alpha_X.$$

## Negative part

In the proposed model, we allow $n_{x_i}$, $i = 1, \dots, 5$, to be positive or negative. If one of these parameters is negative, the associated equation can be written equivalently as, for $t \ge 0$,

$$\frac{dx_i}{dt} = \frac{\tilde{k}_{x_i} E_y^{m_{x_i}}}{b^{m_{x_i}} + E_y^{m_{x_i}} x_i^{m_{x_i}}} x_i - \mu_{x_i} x_i,$$

with

$$m_{x_i} := -n_{x_i} > 0 \quad \text{and} \quad \tilde{k}_{x_i} := k_{x_i} b^{m_{x_i}}.$$

<figure><img src="image_4.png" /><figcaption>Figure 3: The case where $n$ is negative: Hill function $y \mapsto f_x(y) = \frac{ky^n}{b^n+x^n y^n}$ with $b = 2$ and $k = 1$. In this case, the function $f_x(y)$ is decreasing with respect to $y$. Left: influence of $n$ (negative) with $x = 1$. Right: influence of the element $x$ with $n = -1 < 0$.</figcaption></figure>

## Stiffness-dependant model

So far, the variations in protein concentration and mineralization depend on time. It is crucial to emphasize at this juncture that the utilization of time-varying differential or partial differential equations characterizes most models about the evolution of proteins over time. These equations are subsequently juxtaposed with experimental time-series data to elucidate and validate parameter values. In our specific circumstance, this methodology was not applicable due to the restriction of having access to data at only a single time point. Consequently, employing time-varying differential equations appeared incongruous. Nevertheless, as variations in protein concentration were provided as a function of changes in stiffness as mentioned in the Experimental Result subsection and illustrated in Fig. 4, we were compelled to devise a means of transmuting our time-varying equations into equations contingent on stiffness. The ensuing sections expound precisely upon this method. To obtain an equation showing variations with respect to stiffness, we proceeded in the following way: to get rid of the time variation from the System (2), we divided each side of the equation in
$x_i$ by the equations in $E_y$. For technical reasons, we assume that the output of two equations is bounded and that the output of the stiffness equation is different from zero. We end up with only one differential equation of $x_i$ with respect to the stiffness given by

$$ \frac{dx_i}{dE_y} = \frac{(k_{x_i} x_i^{n_{x_i}+1} - \mu_{x_i} x_i (b_{x_i}^{n_{x_i}} + x_i^{n_{x_i}} E_y^{n_{x_i}})) (b_S^{n_S} + \sum_{i=1}^5 (\alpha_i x_i^{n_S} E_y^{n_S}))}{(k_S E_y^{n_S+1} - \mu_S E_y (b_S^{n_S} + \sum_{i=1}^5 (\alpha_i x_i^{n_S} E_y^{n_S}))) (b_{x_i}^{n_{x_i}} + x_i^{n_{x_i}} E_y^{n_{x_i}})} \quad (3) $$

The initial condition for stiffness is equal to 8 kPa for each GRN element $x_i$ (see Fig. 4). System (2) is a system of two equations depending on $i$, which we consistently refer to as System (2). For Equation (3), it actually represents several equations depending on $i$, but we refer to it in the text as Equation (3). It is formally assumed that the denominator in Equation (3) is nonzero (see the second point of the Remark below).

## Remarks.

1. Our decision to treat stiffness not as a fixed quantity, but as an evolving one (hence the use of the variable $y$ in our equations), is motivated by the fact that stiffness changes over time due to the progressive mineralization of the collagen matrix. As osteoblasts release varying concentrations of proteins $x_i$, the matrix stiffens, making stiffness a time-dependent quantity linked to protein dynamics. However, since experimental data only provide a single time point ($t = 336$ hours) without initial conditions, fitting time-dependent simulations to Figure 4 is not feasible. We therefore adopted a reverse approach: each experimentally measured stiffness value is assumed to elicit a specific osteoblastic response via protein release. This assumes a closed (autonomous) system with a feedback loop, where osteoblast mechanosensors detect and respond to stiffness. Our model thus emphasizes bidirectional interactions—stiffness influences osteoblasts, which in turn modify stiffness—rather than proposing a purely phenomenological, time-driven description.

2. Another possibility would be to include additional events or processes in this complex mechanism—something that could be biologically relevant. However, our study focuses on a closed system as a first approach. The inclusion of external factors such as time-dependent events, pathological conditions, or other regulatory mechanisms will be addressed in future work.

3. We note here that System (2) as well as the fact that $E_y^* = 0$ is one of the steady states and using the Cauchy-Lipschitz existence and uniqueness theorem, it is quite straightforward to justify that the first equation of System (2) (which appears in the part of the denominator of Equation (3) - showing a negative term) cannot be equal to zero when initial conditions are positive (which is our biological assumption). Furthermore, it is straightforward to show that the solution remains positive. On the other hand, if a positive equilibrium exists, it could cause the denominator in question to vanish. Given the highly nonlinear nature of the system, it is difficult to analyze this issue analytically, and we leave it for future work. However, numerical results suggest that everything behaves correctly within the interval considered: the derivative of $x_i$ with respect to $E_y$ remains well-defined and bounded.

The next step is to approximate the solution by an ODE scheme and then used the Matlab® command `lsqcurvefit` to estimate the parameters.

## Results

The results section of this article presents the findings of our investigation into the impact of microenvironment stiffness on bone cell behaviour and protein production. Through a combination of in-vitro experiments culturing bone cells on different stiffness substrates, adapting a previously
published model from being continuously time-dependent to continuously stiffness-sensitive, and
simulating protein production in relation to stiffness and other protein production based on the best
estimate of parameters obtained from the experimental work, we were able to collect the following
Experimental Results and Simulation Results.

## Experimental Results

Our experimental findings, as illustrated in Fig. 4, delineate a distinct pattern in protein production
in response to varying microenvironmental stiffness. Specifically, we observed an augmentation in
the production of RUNX2, OSX, and OC with increasing stiffness, while BSP and OPN exhibited a
decrease in production under the same conditions. These stiffness-sensitive trends were discerned at
a timepoint of *t* = 336 hours (that is two weeks), affirming the pivotal role of the microenvironment's
physical attributes in regulating bone cell behaviour and protein production. The observations vali-
date that heightened stiffness corresponds to increased RUNX2, OSX, and OC production, peaking
at 64 kPa. Conversely, BSP and OPN showcase reduced production under heightened stiffness, with
BSP reaching its zenith at 16 kPa and OPN at 8 kPa. Examination of RUNX2 reveals maximal
expression at 64 kPa, contrasting with minimal expression at 8 kPa. In the case of OC (Osteo-
calcin), the pinnacle of expression occurs at 64 kPa, while minimal expression is noted at 16 kPa.
Simultaneously, OSX (Osterix) manifests maximal expression at 64 kPa and minimal expression at
16 kPa.

## Statistical analysis

PCR data were normally distributed for BSP, OPN, RUNX2, and OSX and not normally distributed
for OC and stiffness; no parametric relationship was found, non parametric relationship between
RUNX2 and stiffness was found (Spearmanρ = 0.713, p = 0.009).
<figure><img src="image_5.png" /><figcaption>Figure 4: Production of BSP, OPN, RUNX2, OC, and OSX under different stiffness values from 8 to 64 kPa, achieved through various coating methods. The plastic substrate, corresponding to the uncoated petri dish, is also reported; while its stiffness is not experimentally controlled, it typically ranges between 2 and 4 GPa. The results are taken at t = 336 hours (= two weeks). Error bars depict a single standard deviation of the mean data acquired per RNA purification and quantification (relative fold change) for the RNA transcripts of bone-specific proteins and transcription factors analyzed in this study.</figcaption></figure>

Simulation Results

The simulation results in Fig. 5 show a remarkable agreement with the experimental data collected given in Fig. 4. We simulated protein production in response to varying stiffness continuously. For the simulations, we solved the nonlinear ordinary differential equations (3) using an explicit Runge-Kutta method. The state variables $x_i$, $i = 1, 2, 3, 4, 5$, representing elements of the GRN, were approximated with respect to the stiffness $E_y$ considered as a variable here. Parameter estimation was done using the least squares method. This standard method allowed us to have the curve that

270

271

272

273

274

275

276
<table><caption>Table 2: Table of the model parameters used for simulation presented in Fig. 5.</caption><thead><tr><th>Stiffness</th><th>Runx2</th><th>Bsp</th><th>Oc</th><th>Opn</th><th>Osx</th></tr></thead><tbody><tr><td>k<sub>S</sub> = 0.9419</td><td>k<sub>R</sub> = 0.8907914</td><td>k<sub>B</sub> = 0.6963</td><td>k<sub>C</sub> = 1.4582</td><td>k<sub>P</sub> = 2.248899</td><td>k<sub>X</sub> = 0.167799</td></tr><tr><td>μ<sub>S</sub> = 0.052096</td><td>μ<sub>R</sub> = 0.5226149</td><td>μ<sub>B</sub> = 1.139299</td><td>μ<sub>C</sub> = 0.698099</td><td>μ<sub>P</sub> = 0.022901</td><td>μ<sub>X</sub> = 0.0792</td></tr><tr><td>b<sub>S</sub> = 1.000099</td><td>b<sub>R</sub> = 1.0046971</td><td>b<sub>B</sub> = 1.113299</td><td>b<sub>C</sub> = 0.965999</td><td>b<sub>P</sub> = 1.151899</td><td>b<sub>X</sub> = 0.0918</td></tr><tr><td>n<sub>S</sub> = 1.230899</td><td>n<sub>R</sub> = <b>-0.4711</b></td><td>n<sub>B</sub> = 0.6594</td><td>n<sub>C</sub> = 0.0462</td><td>n<sub>P</sub> = <b>-1.2103</b></td><td>n<sub>X</sub> = 0.1706</td></tr><tr><td></td><td>α<sub>R</sub> = 1.000599</td><td>α<sub>B</sub> = 1.012499</td><td>α<sub>C</sub> = 0.996599</td><td>α<sub>P</sub> = 1.052399</td><td>α<sub>X</sub> = 0.099299</td></tr></tbody></table>

as close as the data by minimizing the sum of the squares of the residuals. Here, the residuals are the differences between the observed values (true data points) and the predicted values produced by the model. Our findings were consistent with the trends observed in the *in-vitro* experiments. Specifically, the model accurately predicted increased RUNX2, OSX, and OC production and decreased BSP and OPN as the microenvironment stiffness increased. These results suggest that our model effectively captures the underlying mechanisms that govern the behaviour of bone cells in response to mechanical stimuli. Regarding the resulting simulations with respect to stiffness, our results indicate that OPN reaches its maximum at a low stiffness of 8 kPa. At the same time, RUNX2, OSX, and OC exhibit their maximum production at the highest stiffness value of 64 kPa. Interestingly, the maximum production level of BSP occurs at an intermediate stiffness value of approximately 11.3 kPa, which presents a unique and intriguing pattern compared to the other proteins (Table 3). In addition, we conducted a 3D phase space simulating 2 protein interactions with respect to stiffness (Fig. 6) and focused on the effects of varying stiffness and protein production levels on bone cell behaviour. By utilizing our stiffness-sensitive model and experimental data, we predicted the production of various proteins and transcription factors, including RUNX2, OSX, OC, BSP, and OPN, in response to changes in both stiffness and other protein productions. We unfolded one example of OPN-BSP-OC interactions (Fig. 7).

In the context of our study, it is crucial to clarify that our approach does not align with a typical validation process. Instead, the chosen parameter estimation method was meticulously designed to minimize differences between computational simulations and real-world observations. Consequently, the agreement observed when comparing experimental results with the model's predictions was expected. To enhance the study's rigour, two main strategies were employed. Firstly, a new data point was introduced to strengthen the validation, with further details provided in the recommended section and subsequent discussion. Secondly, a precise parameter fitting approach was implemented, calibrating all stiffness coefficients except for a single instance. Looking ahead, future research should consider expanding the dataset beyond conventional benchmarks and exploring alternative modelling methods to deepen the findings.
<figure><img src="image_6.png" /><figcaption>Figure 5: The resulting simulations for stiffness and estimation of parameters. The curves are drawn from Equation (3). These graphs visualize the RNA transcripts of protein levels at different stiffness values. The experimental data are represented as dots with the error bars that depict a single standard deviation of the mean data acquired per RNA purification and quantification (relative fold change).</figcaption></figure>

<table><caption>Table 3: The resulting simulations with respect to stiffness for estimation of parameters, maximum values are in bold. Each row corresponds to a different specific stiffness value ($E_y$), expressed in kilopascals (kPa), while each column represents a specific factor or protein, namely bone Sialoprotein (BSP), Runx2, Osteocalcin (OC), Osteopontin (OPN), and OSX. Values used are in Table 2 with the max values in bold. The production of BSP, Runx2, OC, OPN, and OSX is influenced differently as stiffness changes with certain proteins showing maximum production at higher stiffness levels and others at lower levels.</caption><thead><tr><th>Stiffness (kPa)</th><th>BSP</th><th>Runx2</th><th>OC</th><th>OPN</th><th>OSX</th></tr></thead><tbody><tr><td>8</td><td>6.7856</td><td>1.9464</td><td>3.0188</td><td><strong>48.6461</strong></td><td>2.3479</td></tr><tr><td>11.3</td><td><strong>7.2128</strong></td><td>2.1006</td><td>2.8009</td><td>48.1026</td><td>2.2963</td></tr><tr><td>64</td><td>5.2542</td><td><strong>2.6493</strong></td><td><strong>3.8227</strong></td><td>42.5320</td><td><strong>2.9955</strong></td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>Figure 6: Example of production dynamics of two proteins according to stiffness, red shows the maximum production, and blue shows the minimum production. Each figure represents a different set of proteins. Our findings show the relationships between the production dynamics of all observed elements and the stiffness parameter denoted as Ey, expressed in kiloPascals (kPa), which pertains to the mechanical rigidity of the cell's environment. For example, when examining the third row, it becomes apparent that under similar bone Sialoprotein (BSP) levels and stiffness, Osteocalcin (OC) and Osteopontin (OPN) exhibit contrasting production behaviors. The curves are drawn from Equation (3).</figcaption></figure>

## Remarks.

These simulations provide insight on four key fronts.

1. First, some results reinforce previously reported experimental observations, lending credibility to the model's structure and assumptions.

2. Second, the model enables predictive analysis by extrapolating beyond available data, offering testable scenarios that extend current biological understanding.

3. Third, these predictions give rise to new hypotheses regarding the regulatory mechanisms at play, which could guide future lines of investigation.

304

305

306

307

308

309

310

311
<figure><img src="image_8.png" /><figcaption>Figure 7: An illustration of the unfolding 3D phase plane production dynamics reveals that the production of Runx2 is associated with a decrease in Osteopontin (OPN). In contrast, the production of bone Sialoprotein (BSP) is linked to a decrease in Osteocalcin (OC). The integral curves are drawn from Equation (3).</figcaption></figure>

4. Fourth — and perhaps most importantly — this framework identifies potential avenues for experimental validation, thereby linking theoretical modeling and empirical research in a productive feedback loop.

It is important to emphasize that the model presented here is dynamic and intentionally designed to remain flexible. As new experimental data become available, the model can be readily updated and recalibrated, improving its accuracy while reducing reliance on conjecture. Once appropriately validated, this framework could serve as a versatile tool for exploring a range of biological conditions, including pathological dysregulations, through targeted adjustments of its para-meters. We also acknowledge the limitations of the current version of the model. As a first attempt to capture the complexity of the system under study, it necessarily simplifies some aspects of the biological reality. Nonetheless, it offers a rigorous foundation upon which more detailed and nuanced representations can be built. In this sense, the present work lays the groundwork for a progressive refinement of theoretical models capable of describing biological phenomena that are currently difficult—or even impossible—to observe directly.

## Discussion

Our theoretical approach depicts the mechanism of epistatic genetic interactions that control mineralization through mechanotransduction. Here, we captured the behaviour of the bone GRN portion that contains mineralization according to stiffness. We showed that the epistatic genetic interactions follow a system of nonlinear differential equations. First, by performing *in-vitro* experiments culturing bone cells on different stiffness, we provided data regarding BSP, OPN, RUNX2, OC, and OSX production at $t = 336$ hours (= two weeks); second, we adapted our previously published model - from continuously time-dependent - to continuously stiffness-sensitive; third, we optimized our model based on the novel experimental results and depicted how the production of transcription factors RUNX2 and OSX are stimulated by mechanotransduction and how the outputs of BSP, OPN, and OC are altered.

In our experimental results, we present a novel dataset capturing the key elements of the bone GRN in response to stiffness stimulation of human mesenchymal stem cells. This approach mimics mechanotransduction by using substrates of varying stiffness to activate the WNT-β catenin pathway. However, the statistical analysis provided here and the empirical observation associated do not fully capture the complex production patterns regulated by the bone GRN under stiffness
activation. Indeed, among the tested elements, only RUNX2 demonstrated a qualitative correlation (Spearman correlation) with substrate stiffness. To enhance empirical observation, in our previous work [8], we proposed a condensed, simple, and theoretical version of the bone GRN involved in the mineralization of the bone matrix, and we built a model consisting of a system of nonlinear differential equations. While we depicted the direct interactions between the genes coding for the transcription factors and those coding for the bone proteins (e.g. enhancers or inhibitors of mineralization), we did not take into account the impact of stiffness that is key to mechanotransduction and bone mineralization pattern.

Here, in the initial phase of our investigation, the focus centred on a singular time point. Our validation methodology was crafted with consideration of stiffness dependence, formulated to mitigate discrepancies between computational simulations and empirical observations. The anticipated congruence between experimental outcomes and model predictions was attained through a standard parameter estimation. A novel data point was introduced to augment the study's robustness, and a precision-driven parameter fitting approach was implemented, excluding a sole coefficient. In contrast to antecedent investigations, our parameterization was not merely transposed but underwent refinement based on novel data. We adjusted our parameters using a time series data set and modelled transcription factors and bone protein production through the Michaelis-Menten and Hill functions. Indeed, the general system - to describe the interactions between the stiffness and each element of GRN - is given by a system that uses the Hill function. Compared to arctangent or polynomial modelling where only phenomenological explanation has been established, the Hill function - and its particular case of Michaelis Menten - is one of the most used models in theoretical biology because of its mechanistic understanding [34]. The classic Hill Equation is often used to describe the evolution in time of the biological interactions that exhibit cooperativity among two binding molecules (like the one we used in [8]). However, it appeared here not to be adapted at all to our new experimental data where only a one-time step was measured ($t = 336$ hours) for four stiffness inputs that are 8, 16, 32, and 64 kPa. Indeed, with such a data set, it is impossible to consider our previous approach. To answer this major concern, we adapted our model to stiffness, and we showed the results of simulations of RUNX2, BSP, OPN, and OC. As previously stated, only a non parametric relationship between RUNX2 and stiffness was found (Spearman= 0.713, p= 0.009), experimentally. Regarding the numerical results, the model parameters were estimated using the nonlinear least squares method. The quality of the model fit was assessed visually, with the model demonstrating a strong agreement with the experimental data - staying within the standard deviation. Parameter estimates and model behavior were interpreted based on their ability to reproduce the observed trends. Regarding the maximum production, our results from Fig. 5 depict pretty well that OPN is at its maximum production when the stiffness is low (8 kPa), and the RUNX2, OSX, and OC are at their maximum when stiffness is at their maximum as well (64 kPa). Interestingly, BSP has the most intriguing behaviour since its maximum arrives at around 11.3 kPa. Our results show that at the start of the WNT-$\beta$ catenin activation, we observe an OPN and BSP production enhancer action. For instance, the best candidate could be Micro-RNA 214 [35] which has been shown to act on BSP directly. Regarding other candidates, miR-9 upregulates BSP [36] with Col 1 - which could work with the increase in stiffness - and OC, which contradicts our results. Regarding OPN, which reaches a maximum at low stiffness, MicroRNA-21 is a candidate for OPN production enhancer [37]. It also enhances RUNX2 and OSX, reaching their maximum of 64 kPa. In addition, the Slug gene has been described as positively correlated with osteoblast markers, including RUNX2, OPN, OC, WNT-$\beta$ catenin signalling mediators, and mineral deposition. However, how Slug gene can regulate its sequence through micro-RNA production remains an open question. Another common mechanism to enhance osteoblast differentiation is via miRNA-mediated suppression of Wnt inhibitors (i.e., DKK1 or GSK3$\beta$) or via targeting negative regulators of RUNX2 (i.e., SMAD6 or SMURF2) [38]. Thus, knowing the cascade of positive micro-RNA action remains an open question in osteoblast differentiation.
Furthermore, when we plot OPN production versus RUNX2 production and OC production versus BSP production in Fig. 7, we show that OPN production decreases faster when RUNX2 production gets larger. BSP production increases when OC production decreases. However, when BSP production is too large, OC production no longer decreases. In addition, OC production seems to impact BSP production with a negative feedback-like. We have quantified how stiffness starts the production of the first bone GRN transcription factor (RUNX2) that turns up the production of OSX, which activates SATB2. Then, SATB2 starts the output of the first enhancers of mineralization (BSP) which induces both the second enhancers (ALP) and the two inhibitors (OC and OPN) of bone mineralization. When considering three elements (in Fig. 6), stiffness alterations directly impact OC and BSP production, while it does not seem to impact OPN production drastically. In Fig. 6, we describe the impact of the stiffness on each bone GRN element (positive in green and negative in orange). While it does not have a biological explanation, it can be understood by looking at how we model the canonical mechanotransduction pathway (e.g. WNT-β catenin pathway activated by stiffness). This pathway depicts how bone cell mechanoreceptor stimulation initiates bone mineralization [10, 11, 12, 13, 14]. In our modelling strategy, we globalize the entire pathway in terms of Stiffness, which means that a decrease in stiffness has to be understood in terms of inhibition of the canonical pathway itself. Several molecules have been depicted has WNT-β catenin inhibitors, such as Dkk proteins and sclerostin [39] at the extracellular level, or micro RNAs [40], complicating the selection of one molecule over another.

Our study's limitations fall into three primary categories. First, temporal constraints arise from relying on data collected at a single time point, limiting our ability to capture the dynamic temporal evolution of biological processes comprehensively. Here, we focused on the outcomes after the differentiation steps, rather than during the differentiation phase. To minimize the impact of the differentiation process and ensure that we were observing only the effects of stiffness stimulation, we extended our investigation beyond the previous one-week period. Second, dataset limitations encompass the potential lack of representativeness, urging the need for broader data acquisition initiatives. For instance, we could investigate experimentally Stab2 and ALP to capture the full picture of the bone GRN. In addition, the following genes, Smad6, Sox2, and ATF4, hold significant roles in diverse signalling cascades and cellular functions despite their exclusion from our study focused on the canonical WNT beta-catenin pathway [23]. Smad6's modulation of BMP signalling, Sox2's involvement in embryonic development and pluripotency, and ATF4's contributions to stress responses and metabolic regulation underscore their versatile roles [24, 25]. While our investigation concentrated on the canonical pathway within osteoblast differentiation, these genes offer intriguing avenues for future research to unravel their broader impact on cellular processes. Additionally, the absence of a comparative analysis with alternative methodologies, along with the inherent assumptions and simplifications in our model, indicates areas for future research refinement and exploration. These identified limitations collectively guide our understanding of the study's constraints and underscore areas for ongoing refinement and advancement in the exploration of intricate biological dynamics. In our previous paper [8], we provided the first theoretical evidence of a necessary genetic inhibition for bone mineralization. Compared to empirical evidence, we had to use an inhibition factor in each equation modelling each element of the bone GRN. It revealed the negative indirect interactions from negative feedback loops or the recently depicted micro-RNAs. Here, our new approach concludes that, when considered under the stiffness perspective, and not only on time, the fitting process brings more information and allows us to explore the influence of the time and stiffness on the protein concentrations. When considering that stiffness comes from the mineralization process of the collagen matrix, it is possible to simulate protein concentration alterations with respect to time and stiffness, considering continuous feedback control. Thus, we can predict the maximum concentration depending on the stiffness of the matrix, and we can predict the other protein productions - or the extracellular matrix stiffness - while inputting only one data. To summarise our computational findings, they suggest a complex interplay between various elements in the bone gene
regulatory network (GRN) influenced by stiffness alterations. The production of key transcription factors (e.g. RUNX2) is initiated by stiffness, leading to a cascade of events involving other elements such as OSX and SATB2. This cascade ultimately affects the production of enhancers of mineralization (BSP and ALP) as well as inhibitors (OC and OPN). The relationship between OC and BSP production appears to involve a negative feedback-like mechanism. Stiffness alterations directly impact OC and BSP production, with a less pronounced effect on OPN production.

Despite the promising findings, the study's limitations should be acknowledged. The reliance on a single time point restricts the ability to capture the dynamic, time-dependent nature of gene regulation and mineralization. Future studies incorporating longitudinal data will be crucial for validating and refining the model. Additionally, the exclusion of certain key genes and alternative signaling pathways, such as the non-canonical WNT pathway and other miRNAs, points to areas for further exploration. Here, while focusing on osteoblast differentiation from mechanotransduction, there are the broader roles of Smad6 (SMAD family member 6), Sox2 (sex-determining region Y)-box 2), and ATF4 (Activating Transcription Factor 4) in diverse pathways [23, 24, 25]. For instance, Smad6, crucial for BMP signalling and bone morphogenesis, exhibits broader relevance in immune responses. Sox2, pivotal in embryonic development, extends its influence on neural development, tissue regeneration, and cancer progression. ATF4, integral to osteoblast function, plays a significant role in stress responses and metabolic regulation. Thus, expanding the dataset to include more variables and a broader range of stiffness conditions could refine our understanding of bone mineralization and its regulation.

This study holds significant implications for advancing bone research, particularly in understanding the role of stiffness in bone mineralization and gene regulation. By incorporating stiffness as a critical factor influencing the bone gene regulatory network (GRN), our model provides a deeper understanding of the complex mechanisms governing bone formation. The predictive capabilities of our approach are particularly valuable for clinical applications, especially in the development of targeted treatments for bone diseases like osteoporosis, osteogenesis imperfecta, and bone fractures. By simulating how stiffness affects protein concentrations and gene activity, this model could help identify new biomarkers and therapeutic targets, leading to more effective strategies for promoting bone regeneration and healing. The predictive power of our model also has broad applications in the design of biomaterials and implants, particularly in the field of tissue engineering. By accounting for the mechanical environment of bone, we can better design scaffolds and materials that mimic the natural stiffness of bone tissue, enhancing the success of bone grafts and improving patient outcomes. Moreover, the identification of key transcription factors and their interactions with stiffness provides an opportunity to develop therapies aimed at modulating the GRN, potentially improving bone formation and mineralization. Our study also contributes to the growing field of mechanobiology, highlighting the essential role of mechanical cues in regulating cellular behavior. This approach could be extended to other tissues and organs that are sensitive to mechanical forces, such as cartilage, muscle, and vascular tissues, broadening the scope of the study's impact. The integration of micro-RNAs into our model further opens up new avenues for molecular-based therapies. By targeting specific micro-RNAs that modulate gene expression in response to stiffness, we can potentially fine-tune the osteoblast differentiation process and improve tissue healing.

In conclusion, while our study offers valuable insights into the role of stiffness in bone gene regulation, further refinement of the model and experimental validation will be necessary to fully realize its potential. However, the study lays a solid foundation for the development of predictive models that can guide clinical interventions, enhance bone health, and advance tissue engineering technologies.
## Conclusion

Our theoretical framework manages to unravel the complexities of mechanotransduction, providing unique insights into the genetic control of mineralization. Acknowledging limitations, such as reliance on a single time point and dataset constraints, our study highlights stiffness-sensitive patterns in key protein production. While OPN peaks at lower stiffness (8 kPa), RUNX2, OSX, and OC reach maximum levels at higher stiffness (64 kPa). BSP's behaviour adds intrigue, peaking at 11.3 kPa. Despite limitations, our findings open avenues for understanding genetic impacts on bone mineralization and predicting mechanical properties, marking a significant stride in precise regulation at the genetic level. The relationship between OC and BSP production appears to involve a negative feedback-like mechanism. Stiffness alterations directly impact OC and BSP production, with a less pronounced effect on OPN production.

487

488

489

490

491

492

493

494

495

496

497

## Acknowledgments

This work was performed within the framework of the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). This work was partially supported by a grant from the IMU-CDC and Simons Foundation.

498

499

500

501

502

## Credit authorship contribution statement

Authors' contributions: J.-Ph. B. and L. P.-M. conceived of the presented idea, A.C. and L. P.-M. developed the theory and performed the computations, J.-Ph. B. and L. P.-M. verified the analytical methods, K. Y.-H. Carried out the biological experiment, J.-Ph. B. wrote the manuscript with support from L. P.-M., K. Y.-H., and A.C., J.-Ph. B. and L. P.-M. supervised the project. All authors discussed the results and contributed to the final manuscript.

503

504

505

506

507

508

## References

[1] Buenzli P.R. Osteocytes as a record of bone formation dynamics: A mathematical model of osteocyte generation in bone matrix. Journal of Theoretical Biology. 2015;364:418 – 427. doi:doi.org/10.1016/j.jtbi.2014.09.028.

[2] Isaksson H., van Donkelaar C., Huiskes R., Ito K. A mechano-regulatory bone-healing model incorporating cell-phenotype specific activity. Journal of Theoretical Biology. 2008;252(2):230 – 246. doi:https://doi.org/10.1016/j.jtbi.2008.01.030.

[3] Rieger R., Hambli R., Jennane R. Modeling of biological doses and mechanical effects on bone transduction. Journal of Theoretical Biology. 2011;274(1):36 – 42. doi:https://doi.org/10.1016/j.jtbi.2011.01.003.

[4] Fritton S.P., Weinbaum S. Fluid and solute transport in bone: flow-induced mechanotransduction. Annual review of fluid mechanics. 2009;41:347–374.

[5] Fisher S., Franz-Odendaal T. Evolution of the bone gene regulatory network. Current Opinion in Genetics and Development. 2012;22(4):390–397. doi:10.1016/j.gde.2012.04.007.

[6] Hojo H., McMahon A.P., Ohba S. An emerging regulatory landscape for skeletal development. Trends in Genetics. 2016;32(12):774–787. doi:10.1016/j.tig.2016.10.001.
[7] Xu F., Li W., Yang X., Na L., Chen L., Liu G. The roles of epigenetics regulation in bone metabolism and osteoporosis. Frontiers in Cell and Developmental Biology. 2021;8:619301.
doi:10.3389/fcell.2020.619301.

[8] Chekroun A., Pujo-Menjouet L., Falcoz S., Tsuen K., Yueh-Hsun Yang K., Berteau J.P. Theoretical evidence of osteoblast self-inhibition after activation of the genetic regulatory network controlling mineralization. Journal of Theoretical Biology. 2022;537:111005.
doi:https://doi.org/10.1016/j.jtbi.2022.111005.

[9] Komarova K.V., Smith R.J., Dixon S.J., Sims S.M., L.M. Wahl L.M. Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone remodelling. bone, vol. 33, 2, 2206-2015, 2003, doi:https://doi.org/10.1016/j.jtbi.2022.111005.

[10] Chekroun A., Pujo-Menjouet L., Berteau J.P. A novel multiscale mathematical model for building bone substitute materials for children. Materials. 2018;11(6).

[11] Komori T. Signaling networks in RUNX2-dependent bone development. Journal of cellular biochemistry. 2011;3:750-5. doi:10.1002/jcb.22994.

[12] Mullen C.A., Haugh M.G., Schaffler M.B., Majeska R.J., McNamara L.M. Osteocyte differentiation is regulated by extracellular matrix stiffness and intercellular separation. Journal of the Mechanical Behavior of Biomedical Materials. 2013;28:183-194. doi:10.1016/j.jmbbm.2013.06.013.

[13] Robling A.G., Turner C.H. Mechanical signaling for bone modeling and remodeling. Crit Rev Eukaryot Gene Expr. 2009;19(4):319-338. doi:10.1615/critreveukargeneexpr.v19.i4.50.

[14] Yavropoulou M.P., Yovos J.G. Mechanical signaling for bone modeling and remodeling. Journal of musculoskeletal and neuronal interactions. 2016;19(4):221-236. doi:27609037.

[15] Astudillo P. Extracellular matrix stiffness and Wnt/β-catenin signaling in physiology and disease. Biochemical Society transactions. 2015;48(3):1187-1198.

[16] Wittkowske C., Reilly G.C., Lacroix D., Perrault C.M. In vitro bone cell models: impact of fluid shear stress on bone formation. Frontiers in bioengineering and biotechnology. 2016;4:87.

[17] Baron C., Nguyen V.H., Naïli S., Guivier-Curien C. Interaction of ultrasound waves with bone remodelling: a multiscale computational study. Biomechanics and Modeling in Mechanobiology. 2020;19:1755-1764.

[18] Herrmann M., Engelke K., Ebert R., Müller-Deubert S., Rudert M., Ziouti F., et al. Interactions between muscle and bone—where physics meets biology. Biomolecules. 2020;10(3):432.

[19] Ruggiu A., Cancedda R. bone mechanobiology, gravity and tissue engineering: effects and insights. Journal of tissue engineering and regenerative medicine. 2015;9(12):1339-1351.

[20] Boskey A.L. Noncollagenous matrix proteins and their role in mineralization. bone and Mineral. 1989;6(2):111-123. doi:10.1016/0169-6009(89)90044-5.

[21] Morgan S., Poundarik A.A., Vashishth D. Do non-collagenous proteins affect skeletal mechanical properties? Calcified Tissue International. 2015;97(3):281-291. doi:10.1007/s00223-015-0016-3.

[22] Schweighofer N., Aigelsreiter A., Trummer O., Graf-Rechberger M., Hacker N., Kniepeiss D., et al. Direct comparison of regulators of calcification between bone and vessels in humans. bone. 2016;88:31-38. doi:10.1016/j.bone.2016.04.016.
[23] Wang R.N., Green J., Wang Z., Deng Y., Qiao M., Peabody M., et al. bone Morphogenetic Protein (BMP) signaling in development and human diseases. *Genes Dis.* 2014;1(1):87–105.

[24] Basu-Roy U., Ambrosetti D., Favaro R., Nicolis S.K., Mansukhani A., Basilico C. The transcription factor Sox2 is required for osteoblast self-renewal. *Cell Death Differ.* 2010;7(8):1345–53.

[25] Makowski A.J., Uppuganti S., Wadeer S.A., Whitehead B.J., Rowland J.M., Granke M., Mahadevan-Jansen A., et al. The loss of activating transcription factor 4 (ATF4) reduces bone toughness and fracture toughness. *bone*. 2014;62:1–9.

[26] Sun M., Chi G., Li P., Lv S., Xu J., Xu Z., et al. Effects of Matrix Stiffness on the Morphology, Adhesion, Proliferation and Osteogenic Differentiation of Mesenchymal Stem Cells. *International journal of medical sciences*. 2018;15(3):257–268. doi:10.7150/ijms.21620.

[27] Bouleftour W., Bouet G., Granito R.N., Thomas M., Linossier M.T., Vanden-Bossche A., et al. Blocking the expression of both bone sialoprotein (BSP) and osteopontin (OPN) impairs the anabolic action of PTH in mouse calvaria bone. *J Cell Physiology*. 2015;230(3):568 – 577.

[28] Ikegame M., Ejiri S., Okamura H. Expression of Non-collagenous bone Matrix Proteins in Osteoblasts Stimulated by Mechanical Stretching in the Cranial Suture of Neonatal Mice. *Journal of Histochemistry & Cytochemistry*. 2019;67(2):107–116. doi:10.1369/0022155418793588.

[29] Tu Q., Zhang J., Paz J., Wade K., Yang P., Chen J. Haploinsufficiency of Runx2 results in bone formation decrease and different BSP expression pattern changes in two transgenic mouse models. *Journal of Cellular Physiology*. 2008;217(1):40–47. doi:10.1002/jcp.21472.

[30] Depalle B., Duarte A.G., Fiedler I.A., Pujo-Menjouet L., Buehler M.J., Berteau J.P. The different distribution of enzymatic collagen cross-links found in adult and children bone result in different mechanical behavior of collagen. *bone*. 2018;110:107–114.

[31] Berteau J.P., Gineyts E., Pithioux M., Baron C., Boivin G., Lasaygues P., et al. Ratio between mature and immature enzymatic cross-links correlates with post-yield cortical bone behavior: an insight into greenstick fractures of the child fibula. *bone*. 2015;79:190–195.

[32] Liu Q., Li M., Wang S., Xiao Z., Xiong Y., Wang G. Recent Advances of Osterix Transcription Factor in Osteoblast Differentiation and bone Formation. *Frontiers in Cell and Developmental Biology*. 2020;8:601224. doi:10.3389/fcell.2020.601224.

[33] Ferrell J.E. Tripping the switch fantastic: how a protein kinase cascade can convert graded inputs into switch-like outputs. *Trends Biochem Sci.* 1996;21(12):460–466.

[34] Somvanshi P.R., Venkatesh K.V. In: Dubitzky W, Wolkenhauer O, Cho KH, Yokota H, editors. Hill Equation. New York, NY: Springer New York; 2013. p. 892–895. Available from:

[35] Xin Z., Cai D., Wang J., Ma L., Shen F., Tang C., et al. MiR-214 regulates fracture healing through inhibiting Sox4 and its mechanism. *Journal of Musculoskeletal & Neuronal Interactions*. 2020;20(3):429.

[36] Liu X., Xu H., Kou J., Wang Q., Zheng X., Yu T. MiR-9 promotes osteoblast differentiation of mesenchymal stem cells by inhibiting DKK1 gene expression. *Molecular biology reports*. 2016;43:939–946.

[37] Oka S., Li X., Zhang F., Tewari N., Ma R., Zhong L., et al. MicroRNA-21 facilitates osteoblast activity. *Biochemistry and Biophysics Reports*. 2021;25:100894.
[38] Hensley A.P., McAlinden A. The role of microRNAs in bone development. bone. 2021;143:115760.

[39] Kim J.H., Liu X., Wang J., Chen X., Zhang H., Kim S.H., et al. Wnt signaling in bone formation and its therapeutic potential for bone diseases. Therapeutic advances in musculoskeletal disease. 2013;5(1):13-31.

[40] Lei Y., Chen L., Zhang G., Shan A., Ye C., Liang B., et al. MicroRNAs target the Wnt/β-catenin signaling pathway to regulate epithelial-mesenchymal transition in cancer. Oncology Reports. 2020;44(4):1299-1313.

[41] Sanders R, Mason DJ, Foy CA, Huggett JF Considerations for accurate gene expression measurement by reverse transcription quantitative PCR when analysing clinical samples. Anal Bioanal Chem. 2014 Oct;406(26):6471-83.