Analysis of mechanotransduction dynamics during combined mechanical stimulation and modulation of the extracellular-regulated kinase cascade uncovers hidden information within the signalling noise

Osteoporosis is a bone disease characterized by brittle bone and increased fracture incidence. With ageing societies worldwide, the disease presents a high burden on health systems. Furthermore, there are limited treatments for osteoporosis with just two anabolic pharmacological agents approved by the US Food and Drug Administration. Healthy bones are believed to be maintained via an intricate relationship between dual biochemical and mechanical (bio-mechanical) stimulations. It is widely considered that osteoporosis emerges as a result of disturbances to said relationship. The mechanotransduction process is key to this balance, and disruption of its dynamics in bone cells plays a role in osteoporosis development. Nonetheless, the exact details and mechanisms that drive and secure the health of bones are still elusive at the cellular and molecular scales. This study examined the dual modulation of mechanical stimulation and mechanotransduction activation dynamics in an osteoblast (OB). The aim was to find patterns of mechanotransduction dynamics demonstrating a significant change that can be mapped to alterations in the OB responses, specifically at the level of gene expression and osteogenic markers such as alkaline phosphatase. This was achieved using a three-dimensional hybrid multiscale computational model simulating mechanotransduction in the OB and its interaction with the extracellular matrix, combined with a numerical analytical technique. The model and the analysis method predict that within the noise of mechanotransduction, owing to modulation of the bio-mechanical stimulus and consequent gene expression, there are unique events that provide signatures for a shift in the system's dynamics. Furthermore, the study uncovered molecular interactions that can be potential drug targets.

ED'A, 0000-0003-1471-5077; AS, 0000-0002-6137-8205 Osteoporosis is a bone disease characterized by brittle bone and increased fracture incidence. With ageing societies worldwide, the disease presents a high burden on health systems. Furthermore, there are limited treatments for osteoporosis with just two anabolic pharmacological agents approved by the US Food and Drug Administration. Healthy bones are believed to be maintained via an intricate relationship between dual biochemical and mechanical (bio-mechanical) stimulations. It is widely considered that osteoporosis emerges as a result of disturbances to said relationship. The mechanotransduction process is key to this balance, and disruption of its dynamics in bone cells plays a role in osteoporosis development. Nonetheless, the exact details and mechanisms that drive and secure the health of bones are still elusive at the cellular and molecular scales. This study examined the dual modulation of mechanical stimulation and mechanotransduction activation dynamics in an osteoblast (OB). The aim was to find patterns of mechanotransduction dynamics demonstrating a significant change that can be mapped to alterations in the OB responses, specifically at the level of gene expression and osteogenic markers such as alkaline phosphatase. This was achieved using a three-dimensional hybrid multiscale computational model simulating mechanotransduction in the OB and its interaction with the extracellular matrix, combined with a numerical analytical technique. The model and the analysis method predict that within the noise of mechanotransduction, owing to modulation of the bio-mechanical stimulus and consequent gene expression, there are unique events that provide signatures for a shift in the system's dynamics. Furthermore, the study uncovered molecular interactions that can be potential drug targets.

Background
Osteoporosis as a skeletal disease presents a burden to health systems worldwide while the average life expectancy increases. In the UK, it has been predicted that it will cost the NHS over £2.2 billion per year in 2025, with the expected cost rising annually [1]. The disease causes changes in bone architecture, a reduction in its optimum mechanical properties and thus an increased likelihood of fractures [2]. The disease provides a multiscale challenge as bone mechanical properties and strength, at the tissue level, are a product of intercellular and intracellular interactions. The building unit of bones is the bone remodelling unit (BRU), which consists of osteocytes (OCys), osteoblasts (OBs) osteoclasts (OCs) and the bone extracellular matrix (ECM) [3]. OBs contribute to the anabolic process by depositing organic material such as collagen, alkaline phosphatase (ALP), osteopontin (OPN) and osteocalcin (OCN) in the ECM, leading to the formation of hydroxyapatite and ultimately giving bone its mechanical properties (such as strength and elasticity). They also differentiate into OCys [4,5]. Conversely, OCs mediate the catabolic process by breaking down organic material within the ECM. OCys are widely regarded as the main mechanosensory cell within bone; therefore, they orchestrate the balance between the catabolic and anabolic process via communication with OBs and OCs via ligands such as transforming growth factor beta [3]. It is noteworthy that signalling through sex hormones is believed to play an important role in osteoporosis development, where four out of five cases develop in post-menopausal women and one out of five in ageing men. Hence, OCys are continuously exposed to dual biochemical and mechanical (bio-mechanical) cues. Using intracellular signalling mechanisms, they sense, integrate and respond to these stimuli to maintain the balance between catabolic and anabolic events. Disturbance of this balance leads to osteoporosis.
Integration of the bio-mechanical cues, and their transduction into cellular responses, within the BRU, are conducted by intracellular activities such as the mechanotransduction process (figure 1). Mechanotransduction involves the sensation of mechanical forces through mechanoreceptors such as integrins [6][7][8][9]. Mechanical signals are converted intracellularly into biochemical reactions via pathways such as the extracellularregulated kinase (ERK) [10,11]. The ERK pathway is involved in the expression of ECM proteins that control bone material properties and mineral density (BMD), such as OCN and OPN. Recently, dual bio-mechanical stimulation with intermittent parathyroid hormone (PTH) treatment and mechanical stimulation was shown to increase BMD and bone formation in mice, thus suggesting a promising treatment for osteoporosis [12,13]. Furthermore, it was demonstrated that bone density can be increased before and during the menopause by using mechanical loading regimes such as resistance training, thus limiting the risk of osteoporosis development. However, precise regimens required to induce a potent therapeutic effect are yet to be determined. This is because of an incomplete understanding of the intricate interactions between the different components across scales, and the dynamic stimulations exerted on cells within the BRU [14][15][16][17]. Therefore, there is a need to develop methods and tools to advance our  Figure 1. Simulating dual biochemical and mechanical excitations using a hybrid multiscale model. The model simulates mechanotransduction within a single osteoblast and links the outcome to a mechanical model that handles the biomechanics at the tissue scale. (a) The mechanotransduction cascade was modelled from the integrin mechanoreceptors and involved the recruitment of the Raf-MEK-ERK signalling module in the cytoplasm; this resulted in the activation of nuclear events and the mRNA expression of osteogenic markers such as bone sialoprotein (BSP) and alkaline phosphatase (ALP). mRNA is translated to their corresponding proteins when they interact with activated ribosomes. The mechanotransduction cascade was represented using the System Biology Graphical Notation (SBGN); a detailed diagram is provided in electronic supplementary material, figure S1. (b) A schematic representation of the mechanical model. A mechanical force is generated using a mathematical algorithm that propagates through the tissue composite as shear stress. The mechanical force is transmitted to integrins, causing them to stretch; the calculated force on the integrin determines its activation. The force value is calculated and passed to the agent-based model to determine if integrins become fully activated and thus propagate mechanotransduction. The large arrows indicate the information passed between the two models. The two models are in contact at every iteration, which is equivalent to 1 s.
understanding of the interactions between the different components and consequently assist in the development of novel treatments for osteoporosis [5,18,19]. One avenue towards this is to comprehend the interplay between bio-mechanical cues, mediated through mechanotransduction, and their influence on cellular responses controlling ECM material properties [20,21]. This study presents a computational approach where a hybrid multiscale model, combining a mechanical model with the agent-based approach (Mech-ABM), was used to examine the modulation of bio-mechanical signals and their influence on cellular events leading to the expression of ECM proteins (ECMp) that modify bone material properties (such as ALP and BSP). Simultaneously, an analytical mathematical method was used to extract key information from the noise within mechanotransduction to elucidate the significant impact of the dual alteration on ECMp expression. This study demonstrates that the dual modulation of mechanical stimuli and mechanotransduction dynamics, at the level of the ERK pathway, impact modalities of ECMp mRNA expression and that the activation state of integrins may contribute to the shift of the system dynamics, which ultimately contribute to the emergence of osteoporosis [17,20,21]. Furthermore, Mech-ABM and the analysis method further emphasize the importance of interfering with interactions within the Raf-MEK-ERK (rapidly accelerated fibrosarcoma-MAPK/Erk kinase-extracellular signal-regulated kinase) module to develop candidate pharmaceutical molecules which add potency when combining bio-mechanical stimulations for future osteoporotic therapies.

Methods
The study examined the impact of dual modulation of biomechanical stimulations on molecular and cellular responses which impact ECM composition and ultimately its material properties. These are centred around the modulation of gene expression and translation of ECMp. This was achieved via utilization of a multiscale hybrid mechanical agent-based model (Mech-ABM) (figure 1) as published previously [8]. Briefly, the ABM was used to mimic the dynamic interplay between mechanical stimulation and intercellular response within a three-dimensional spherical OB embedded within a non-demineralized bone ECM. The ABM simulated mechanotransduction downstream of the integrin mechanoreceptors, which recruited the focal adhesion protein (FAK) and the ERK signalling pathway and led to the synthesis and deposition of ECM proteins (figure 1a) [8]. The cell was separated into two main compartments, the cytoplasm and the nucleus, which were separated from each other and the ECM by membranes. The agent molecules, each of which is represented as a discrete interacting entity, were assumed to be homogeneously distributed inside their respective compartments, and their movement was described via Brownian motion. These agents (including the cell, see electronic supplementary material, table S1) were communicating X-machines, which rely on state transitions, memory and transition functions to execute processes. Electronic supplementary material, figure S1 illustrates a representation of all the molecular interactions and types of molecules simulated. The ABM did not include sophisticated organelles such as the Golgi apparatus.
Molecule-molecule interactions, binding events and activation-deactivation state transition (via activation cycle switch (ACS)) are the driving mechanisms of ABM. Furthermore, the Mech-ABM integrated stochastic processes that updated every agent's global and local variable over time; specifically ACS [22]. These were derived from probability functions Ψ(t) [22][23][24] The mechanical model mimicked, at the tissue level, the tissue biomechanics and associated mechanical perturbation, as detailed previously [8]. Briefly, a mechanical shear stress was applied at an infinite three-dimensional composite representing the bone ECM; the mechanical forces were transmitted through the ECM to the ECM-cell interface, where integrins (represented as springs) formed the connection between the ECM and the cell. The mechanical force stretches the integrins and, at a particular mechanical threshold, activates them, hence initiating intercellular mechanotransduction. The accumulated number of protein agents and mRNA agents over time (24 h) partitioned by their state variables (e.g. active/inactive and bound/unbound) was the primary output analysed from the Mech-ABM. The external mechanical perturbation was a periodic square wave characterized by the amplitude M and the oscillation period P. Modulations of mechanical cues were examined by varying the M (M = 100 and 10 000 µPa) and P (P = 1000, 5000 and 20 000 s). Within the range of values of P investigated, two were much shorter than the total duration of the sampling period, named the epoch, and, on one occasion, the period was set so that the stimulation applied was constant for the entire simulation. The biochemical perturbation was achieved via modulating intercellular molecular interactions within the ERK pathway associated with the protein activation cycle (i.e. feedback loops). The total number of simulations was 1800, with 15 repeats per condition.
Re-scaling the three-dimensional interaction volume of agents in the system and the number of agents, accordingly, has little impact on the overall rates of intracellular reactions within an ABM framework, if the system does not become excessively sparse and there is enough statistics as demonstrated by Rhodes et al. [25]. Hence, the size of a large system of molecules, such as an OB, can be reduced without losing the complexity of the dynamics.
From the simulated processes, the resulting signals representing the number of molecules in a given molecular state show a series of spikes ( prominent peaks), that is, rapid abrupt variations in the signal, which appear to be unsynchronized with the external mechanical perturbations. Nor do the times in between these spikes present a narrow distribution, a typical signature of periodicity, and neither are they synchronized by repeating the stimulation periodically multiple times with identical mechanical loads. Nevertheless, such spikes are characterized by amplitudes of several molecules of magnitude, which imply a form of synchronization despite the non-periodicity. In the long run and for an infinite number of molecules, these fluctuations may asymptotically disappear; yet, there exists a scale where they cannot be neglected and cooperative molecular phenomena have various biological functions [14,16]. However, information regarding the complexity of a non-ergodic system and pattern of activation can be obtained. This was achieved by using an analytical method utilizing subordination theory [26,27], as published previously [28]. Briefly, the signal (fluctuating activation spikes) was transformed, via a filter (equations (2.1) and (2.2)), owing to the inherent stochasticity of the system and therefore the ensuing temporal variations in the output (fluctuating signals), whereÂ DT is the moving average operator over a time period [t, t þ DT] and y i (t) is the averaged signal with the reduced level of noise.
royalsocietypublishing.org/journal/rsfs Interface Focus 11: 20190136 A numerical method using subordination theory was used to disentangle random processes into their respective parent process and directing process in order to analyse the latter in terms of patterns of recurrence of events [26]. This was achieved by applying a subordination process separating the natural time and the physical time (the macroscopic and experimentally observable time). A moving average filter with a Hilbert transformation was applied twice, thus determining the signal time series and amplitude as positive definite time series with enhanced prominent peaks and dumped valleys close to the zero axis. The peaks are detected as events ε, if the signal minus one standard deviation σ i of the fluctuations drops after each peak, p n (t), below the corresponding preceding peak value, Therefore, critical events can be detected and provide detail on the repetition of large oscillations within the number of molecules that rapidly reach or depart from a given node within the network. The probability density obtained by normalizing the integral distribution to 1 was expressed as a waiting time distribution (WTD). WTDs between peaks or the time extent of peaks are considered. WTDs are built by sampling the occurrence of peaks for an epoch of time much larger than the average time of two consecutive peaks. A kernel density estimator (KDE) is defined as where w is the support of the operator. Owing to the stochasticity within the model and the applied mechanical stimulation, the critical events would approximately result in a complex overlapping of convolution of exponentially distributed processes characterized by multimodality. WTDs are dependent on the selected parameter for analysis, thus they reflect the specific dynamics of Mech-ABM. The integrated nature built into the Mech-ABM and that of the mechanical excitation result in the critical events being gamma distributed and typified by modalities. This distribution is used to derive signal interim periods (SIPs) τ between consequent events with the magnitude of the fluctuation. The area τ multiplied by the magnitude measures the minimum time cumulated by a particular molecule, at a given activation state, and the length of time it is involved in the same critical event. The total number of molecules within the same critical event over the entire duration of the sampled signal was demonstrated by the non-normalized distribution. Subsequently, the WTD ψ(τ) is the probability density obtained by normalizing the integral distribution to 1. In this study, the WTDs of integrins, ECMps, mRNAs and ECMps were illustrated as they are considered to be primary osteogenic markers in addition to being among the primary components that reflect ECM alteration in response to dual modification in bio-mechanical signalling [4,29,30].

Simulations and sensitivity analysis
Creating the dual modulation of bio-mechanical cues on the dynamics of mechanotransduction and gene expression events were investigated by modifying the external mechanical perturbation and intracellular feedback loops of the ERK pathway. With respect to mechanical excitation, this was via changes to the magnitude (M) and oscillation period (P), while for the feedback loops it was achieved by modulating agents activation cycle (AC) times T • → • (• → • represents the initial and final transition states) for: For repeatability, the model was simulated for 10 independent runs, though with identical initial conditions for each simulation, the distribution was computed and the 10 independent results were used to generate confidence intervals. The top and bottom of the confidence intervals around the WTDs are not probability functions. The generic ABM platform FLAME (flexible large-scale agent modelling environment) was used to simulate intercellular mechanotransduction and mechanical loading [31][32][33][34]. At t 0 , Mech-ABM contained 8890 agents and ran with a unit time step corresponding to 1 s for approximately 24 h. The agents' three-dimensional coordinates within the cell and velocity at t 0 were randomly selected from a uniform distribution. The mechanical load applied at the tissue level was a periodic square wave function defined by its magnitude M and oscillation period P, the minimal magnitude was 100 µPa and the phase was set to zero (see electronic supplementary material,  The SIP range between τ 1 and τ 2 is characterized by multiple peaks changing with the age of the system, as illustrated in figure 4a. Except for T MEKact + ERKd → ERKact < 300 s, in the WTDs of free mRNAs, there is a stable peak around τ 4 = 100 s and an oscillating tail which decays slowly until the next stable mode at τ 3 . Instead, for values of T MEKact + ERKd → ERKact smaller than 300 s, there is no stable peak at τ = 100 s and the decaying tail becomes a long fluctuating plateau (figure 6b). Consequently, faster activation of ERK induces a more focused transcription of mRNAs, while at larger values of T MEKact + ERKd → ERKact there is a wide set of SIPs between τ 1 and τ 2 with high probability that they cause non-rhythmic rapid variations of mRNA production. It is important to emphasize that consistency of behaviours among all types of mRNA, which in the ABM are identically interacting molecules and are dynamically indistinguishable, showing the reliability of the fluctuation analyses and the robustness of the method.

Shifting integrin activation dynamics and their
response to mechanical stimulation is dominated by mechanical stimulation frequency (P) and then by extracellular-regulated kinase activation cycle switch With respect to mechanical perturbations, for all the analysed epochs and the simulated periods (P), the results supported that mechanical perturbation was a key trigger of events,   with integrins as the main mediators. As the mechanical stimuli cease, eventually signal propagation terminates and the number of activated molecules approaches zero. For active integrins, increasing P by a multiple of 40 resulted in the reduction of the variability among independent repetitions of the simulations. In the first and second rows of figure 5, the WTDs with P = 5000 s present a tail with a large confidence interval, and at large τ the shapes largely fluctuate among different epochs. Instead, for P = 200 000 s, the waiting times produce clear bimodal distributions with reduced confidence intervals at large τ, resulting in no tail and a stable shape during different ages. The reduction of noise on active integrins' SIPs is strongly sensitive to the length of the period P, but not to the magnitude of the perturbation (figure 5b , T ERK act → ERK d and T ERK act + RUNX2 d → ERK act larger than the baseline values, the WTDs' noisy tails disappear, and the distributions are the same for any value of the perturbation period.

Ageing characteristics are observed at the level of nuclear events
The system presented previously demonstrated an ageing characteristic [28] at the level of the ERK pathway. This ageing is a feature that is a direct consequence of the initial conditions and corresponds to a slow variation of the WTD as time passes, Viz., where initial conditions and perturbations (such as those of P) might impact the emergence and frequency of critical events monitored, and hence WTDs. Among all   royalsocietypublishing.org/journal/rsfs Interface Focus 11: 20190136 molecules analysed, transcribed factors are characterized by ageing with an evident and clear progression. From the integral distributions in figure 6, we observe that, with the increase of age, there is a reduction in the number of abrupt fluctuations. A similar ageing effect appears in all the mRNAs bound to ribosomes. However, in many cases analysed, the corresponding WTD does not show sufficient discrepancies between different epochs, signifying that dispersion times decrease by the same amount as time passes, and the ageing effect is reduced to a slowing of the system dynamics. As a direct consequence of the fixed number of total ribosomes simulated in the cytoplasm, ribosomes in their free state also present an ageing effect, but their dispersion times increase with the age of the system. For free mRNA, the variations with the age t of the dispersion times change at different rates depending on the SIP value τ, so their WTD does not preserve the shape but changes from a trimodal distribution to a bimodal distribution as time passes (figure 4).

Discussion
The study shows that it is feasible to extract informative detail from the naturally chaotic mechanotransductionmediated gene expression events in OBs. This information can provide additional insight regarding the shift in the system's dynamics owing to modulation of bio-mechanical cues. It further supports that modulation of the molecular interactions of the Raf-MEK-ERK module can provide innovative targets for the development of new drugs for osteoporosis. These outcomes were feasible via utilization of the hybrid multiscale Mech-ABM and an analytical method that used subordination theory.
Osteoporosis develops as a result of a reduction in the efficiency of the anabolic mechanisms that maintain ideal bone material properties and strength. This is directly related to the OBs' response to a bio-mechanical stimulus, where OBs deposit ECMps and differentiate to OCys. Mechanotransduction plays an important role in mediating these events. Hence, inducing OBs towards osteogenic response at the level of mRNA transcription and translation is important. How dual bio-mechanical stimuli achieve this is poorly understood. Nonetheless, we previously demonstrated that fluctuations in the system (OBs' mechanotransduction) are beyond mere system randomness [35], and that, with dual alteration of the mechanical load and ERK's feedback loops (i.e. ACSs), the WTDs can reflect molecular signatures of the system with respect to cytoplasmic events. In this study, increasing ACS caused a shift in the modalities of free mRNA, making gene transcription events more chaotic and frequent, as illustrated with reduced levels of ALP mRNA, increased WTD modalities and an expanded tail. This demonstrates that the modulation of ERK activity can shift the cellular response and enhance the expression of ECMp. This observation agrees with the proposals that Yang et al. [11] recently demonstrated that the ERK cascade is important in orchestrating mechanotransduction and is vital for the maintenance of ECM within bone. Additionally, this impact on osteogenic gene expression is in line with in vitro studies showing that inhibition of ERK activation antagonizes RUNX2 nuclear activity and thus reduces the expression of ALP and OCN, ultimately reducing ECM mineralization [28]. Tomida et al. [36] demonstrated similar behaviour in inflammatory systems, whereby they showed that asynchronous oscillatory activation of the MAPK protein and ERK relative, p38, generated a burst of activity of inflammatory genes. Furthermore, it was demonstrated previously that modulated stochastic and pulsatile activity of ERK could modify homeostatic gene expression events and drive cell proliferation when nudged by external signals [37][38][39]. Consequently, this further suggests that mechanical signals can impose a nudge effect on ERK signalling towards anabolic activity by stabilizing the WTD.
The results presented here illustrate that the ensuing gene expression events of bio-mechanical signals, in contrast with  Figure 6. BSP mRNA bound to ribosome. Spectrogram of the critical events for BSP mRNA during translation. On the x-axis, the age t of the non-normalized histogram also corresponding to the beginning of an epoch; on the y-axis, the SIP τ; on the z-axis, the number of events registered during each epoch of length 10 000 s. The vertical black lines crossing the surface show the 95% CI. The colour defines the value of the WTD ψ at SIP τ and age t. The WTD is bimodal and the shape is constant at all t. Nonetheless, the total number of critical events decreases with age.
royalsocietypublishing.org/journal/rsfs Interface Focus 11: 20190136 cytoplasmic mechanotransduction events [8,35], are noisier. These chaotic events were particularly evident at the level of transcribed ECMp mRNAs (e.g. free ALP mRNA, figure 4) and translated mRNA (figures 2 and 3). The latter was measured by the number of activated mRNA bound to ribosomes. The eruption of activity at the gene expression level, following conservative activation activity in the cytoplasm by the ERK module, is in line with physiological observations [40,41]. Furthermore, our observation that the period of mechanical stimulation causes a shift in the dynamics of the system, by shifting the WTD modality from one value to another, is equivalent to in vitro studies demonstrating that short bursts of mechanical stimulations can elicit a potent cellular response [42,43]. This is probable considering that the most significant impact on mRNA magnitude and WTD modalities was observed with dual manipulation of mechanical load and ERK's ACS. Increasing P and ERK ACS has resulted in a more pronounced shift in the WTDs' modalities, whereby multi-modalities shift to either tri-or bi-modalities. This suggests that dual change of bio-mechanical stimulation drives the system, and in particular ECMp mRNA levels, from somewhat chaotic and frequent to stable and specifically expressed at particular SIPs. This suggests that mechanical stimulation in the form of P acts as a coincidence detector and filter, selecting ERK-induced SIPs that will propagate optimum mRNA transcription and translation events for osteogenic markers such as ECMps, thus substantially augmenting ECMp deposition. These decisions ultimately influence bone density and BMD. This is in line with the work of Albeck et al. [37], who demonstrated that, within genetically identical sister cells, bursts of asynchronous ERK activities determine if cells entered S-phase and thus the dynamics of cell proliferative response to epidermal growth factor stimulation.
As it was demonstrated previously that biased PTH receptor (PTHR) agonists recruit the Raf-MEK-ERK module and enhance its signalling impact [44][45][46], we provide another avenue to explore new therapies where protein interactions within the module can be altered using small peptides. The candidate molecules can be scaffold proteins such as KSR or beta-arrestins. The latter has links with PTHR and was shown to be the mediator of the biased agonism [44,45]. This is a credible suggestion since the use of small molecules aimed at intracellular molecular targets has been explored recently and showed a degree of promise [47][48][49].

Limitations
We acknowledge that our Mech-ABM uses assumptions and constraints that will put some restrictions on what to infer from the data. For example, mechanical stimulation does not change the spherical shape of the cell over the course of the simulation, and this might impose a restriction on the propagation of the mechanotransduction. We also recognize that, to draw a full conclusion on ECM health over time, the model should be simulated for a longer period, at least a few days, to see the long-term effects of dual bio-mechanical stimulation on cell behaviour. However, this is computationally costly and will require adjustment of the model to overcome this issue. Furthermore, we appreciate the role of sex hormones and their involvement in osteoporosis progression; however, these are outside of the scope of the Mech-ABM as its cell signalling network does not include the hormonal component. Yet incorporating these into the model in the future is uncomplicated and is achievable. Though we realize there are limitations, nonetheless, thus far the model has replicated physiological phenomena and the emergent behaviour of bio-mechanical stimulation consistent with physiological observations both in vitro and in vivo. These give us great confidence in the interpretations made using it.

Conclusion
We observed large fluctuations in the information hidden in the noise which are beyond the dynamic variations of molecular baselines. These were observed as alterations in WTDs of free and activated mRNA for osteogenic markers (e.g. ALP and BSP). Therefore, WTDs of each molecule are a signature of the system's dynamics. In contrast with more traditional approaches, where noise has been adopted as a measure to quantify the error around the expected values of time-dependent signals, our study explored and analysed abrupt fluctuations departing from the trend, thus hidden information characteristic of mRNA molecules and of the process's dynamics was extracted. These forecast that protein-protein interactions at the level of the RAF-MEK-ERK module play a key role in dual biomechanical signalling. Thus, the modulation of scaffold proteins could be an innovative approach to modify mechanotransduction and ultimately develop innovative therapies for osteoporosis.
Data accessibility. Our code, data and appropriate materials are available under the Creative Commons Attribution License. These are available on Google Drive with the following hyperlink: https://drive. google.com/drive/folders/0AEwfzIdJmd_BUk9PVA.