ABM of osteoblast's mechanotransduction pathway: time patterns of critical events
Gianluca Ascolani, Timothy M. Skerry, Damien Lacroix, Enrico Dall'Ara, Aban Shuaib
aa r X i v : . [ q - b i o . S C ] F e b Ascolani et al.
ABM of osteoblast’s mechanotransductionpathway: time patterns of critical events
Gianluca Ascolani , Timothy M. Skerry , Damien Lacroix , Enrico Dall’Ara and Aban Shuaib AbstractBackground:
Mechanotransduction in bone cells plays a pivotal role in osteoblast differentiation and boneremodelling. Mechanotransduction provides the link between modulation of the extracellular matrix andintracellular actions. By controlling the balance between the intracellular and extracellular domains, themechanotransduction process determines optimal functionality of the skeletal dynamics, and it is one of thepossible causes of osteophatological diseases.
Results:
Mechanotransduction in a single osteoblast under external mechanical perturbations has beenmodelled in the agent based framework to reproduce the dynamics of the stochastic reaction diffusion processamong molecules in the cytoplasm, nuclear and extracellular domains. The amount of molecules andfluctuations of each molecular class has been analysed in terms of recurrences of critical events. A numericalapproach has been developed to invert subordination processes and to extract the directing processes from themolecular signals in order to derive the distribution of recurrence of events.
Conclusions:
We observed large fluctuations enclose information hidden in the noise which is beyond thedynamic variations of molecular baselines. Studying the system under different parametric conditions andstimuli, the results have shown that the waiting time distributions of each molecule are a signature of thesystem’s dynamics. The behaviours of the molecular waiting times change with the changing of parameterspresenting the same variation of patterns for similar interacting molecules and identifying specific alterationsfor key molecules in the mechanotransduction pathway.
Keywords: osteoblast mechanotransduction; molecular network; fluctuations; subordination theory; directingprocess; waiting time distribution scolani et al.
Page 2 of 31
Background
Mechanotransduction drives the orchestration between bone remodelling, the pro-cess of bone adaptation over time, and the activity of the bone cells: osteoblasts,osteoclasts and osteocytes [1, 2]. Transmission of mechanical stimulations throughextracellular matrix (ECM) medium deep inside the cells affecting the intracellularsignalling, both in the cytoplasm and in the nucleus is a phenomenon which, despitesignificant progress in recent years, has been found more complex than it has beenthought [3–7]. Studying the deviations from the coordinated balance of mechan-otransduction signalling in osteoblasts, which can cause osteodegenerative diseasessuch as osteoporosis, osteopetrosis, and osteoarthritis, is fundamental [8–10].Effects of mechanical loads propagate through the bone structure, stimulatingboth the osteocytes with their dendritic-like processes and the lining cells on thebone apposition surface. The production of osteoapatite and factors involved in theconformation of the ECM is a continuous dynamic process occurring in both casesof fracture repair and tissue deformation. The mechanical stimulations propagatedthrough the tissue are sensed by osteocytes and osteblasts via mechanoreceptorssuch as integrins [11].Integrins play a pivotal role in migration as well as in osteogenic differentiation andosteoblast activation [3, 4, 12, 13]. The extracellular domain of integrins providesanchorage and interaction interface with the ECM, while ensembles of integrins,on the intracellular domain, are molecular sensors accumulating the effects of theexternal focal adhesion forces. One of the main aspects of integrin mediated signalsis the cumulation of structural tension in the protein, and its consequent activationis antagonized by a relaxation process [14, 15]. Thus, transduction can be stronglyprolonged if the external force does not produce adequate deformations.Stretching and compression of non homogeneous tissue locally depends on themechanical properties of the extracellular milieu. Considering the limited diffusionof large non soluble molecules in dense ECM, the local properties of the latter aremainly affected by cell secretion of factors which accumulate in the surroundingvolume. In turn, variation of the ECM stiffness changes how a cell senses the forces.Therefore, the signal conveyed by integrins affects bone tissue at different scales[16–18].Mechanotransduction has been investigated with experimental [8, 19, 20] andcomputational [21, 22] approaches; nevertheless, many aspects of this interaction inbone cells have not yet been explored in large part due to the lack of accurate toolsto study molecular events in live cells. Many of the in vitro setups are limited to 2Ddistributions of cells with reduced density of surrounding ECM in order to facilitatestimulations and observations of the system. Measurements based on traction forcemicroscopy are affected by low spatial resolution and constraints in the orientationof applied forces [23].Molecular tension sensors experiments are useful for observation of downstreamchemical signalling concurrent with fluorescence signals. However, spectrum limi-tations can be challenging when downstream signalling is being explored. Further-more, Foster resonance energy transfer (FRET) experiments often limit the numberof interactions, the range of extensions and forces that can be explored [24]. In silicoMolecular Dynamics models are limited to a simulated time which is in the order of scolani et al.
Page 3 of 31
100 nanoseconds [25], while continuous approaches lack the possibility to considerdiscrete and highly heterogeneous systems that can form local recurrent structuresamong the molecules [26].An agent based approach can overcome some of the above mentioned limitationsof present in vitro and in silico approaches. The advantage of using an agent basedmodel (ABM) consists in generating a complex dynamics based on simple rulesfor the physical interactions that would not be easily obtainable in a continuousrepresentation of the system and would be experimentally difficult to control. ABMsare suitable to mimic 3D cellular systems where the dynamics of the signalling andthe occurrence of events inside a cell can be tracked at the molecular scale [27]Furthermore, the heterogeneity of molecular populations responding to directionaland time modulated stimuli can be considered [28] Additionally, these approachescan be adopted to address the effects of non homogeneity in the plasm and nucleusof the cell far from the condition of well mixed systems assumed in the Gillespie classof algorithms [29]. In contrast to MD and ODEs/PDEs methods, the coarsening ofmolecules to single diffusing interacting agents accounts for intracellular stochasticeffects and yet, run large simulations for periods of time corresponding to one day[26, 30].The aim of this study was to use an ABM model of the osteoblast to study theeffect of external mechanical stimuli to the local molecular events within this cell.
Implementation
We used an ABM to mimic the dynamics in the mechanotransduction pathway ofa single osteoblast during the deposition of ECM factors under external mechanicalperturbations in order to analyse fluctuations of molecular signals and to extracttime patterns of critical events characteristic of the process dynamics. In the ABM,each agent represents a specific molecule, which has a relevant role in the processof ECM formation and is included in the mechanotransduction network (see fig. 1).The agents move following a constrained Brownian motion in one of three geomet-rical compartments – the cell nucleus, the cytoplasm, and the extracellular regionsurrounding the cell (orange borders in fig. 1). When the distance between twoagents is shorter than the interaction range R inter , a set of rules depending on theinternal state of the agents is applied to each agent involved in the interaction. Sim-ilarly, single agents can trigger rules due to their internal clock, state or informationexchanged. The actions resulting from the application of rules triggered by an agentmay involve:1 the changing of its own state,2 the self-destruction,3 the production of information causing:(a) other existing agents to change their states, as well as(b) the creation of new agents, or(c) triggering other actions.In terms of agents, the reactions between two molecules generating a molecularcomplex are two types. The first requires the change of state of one agent andthe disappearance of the other; the second consists in the disappearance of bothagents involved and the appearance of a new agent. Reactions referring to a transfer scolani et al. Page 4 of 31 of a molecular messenger are obtained by the change of the states of the agentsinvolved. Dissociations in two or more molecules of complexes regulated by thedelay and molecular persistence in a given state have two ways. One way is thedisappearance of the original agent and the appearance of as many agents as thenumber of dissociated molecules. The other way consists of creating as many agentsas the number of dissociated molecules minus one and the change of state of theoriginal agent.Among the molecules simulated there are proteins (e.g. integrins), ribosomes,mRNA and vesicles. For the sake of simplicity, other elements like amino acids,organelles, molecules and complexes not included in fig. 1 are neglected. The ruleswe used in the proposed model are based on the affinity of the molecules, and theyare derived by the chemical reactions sketched in the molecular network diagram.The updating scheme adopted [31] requires that, at each iteration, all themolecules are swept and checked for possible rules to be applied. The types of phys-ical interaction driving the model dynamics are two: 1) binary binding of moleculesbased on one molecule querying all the affine available molecules in the sphere withradius R inter and then requesting the interaction with the nearest one; 2) degrada-tion, relaxation and unstable complex dissociation activated at the expiration of aninternal timer set to a random time tossed from a probability function Ψ( t ), whenthe complex changes state and decreased of 1 unit of time at each iteration.In the simulated ABM, among various signals, we generated the cumulated num-ber of agents partitioned by their internal biological state variables potentially repre-senting experimental distinguishable molecular subpopulations (e.g. active/inactiveand bound/unbound).Many of the simulated molecules are present, at least, in two complementary iso-forms: active and inactive. In macroscopic steady state conditions, when one isoformproduces a smooth signal close to zero axis interrupted by positive spikes, the othergenerates a negative fluctuation around a positive nominative value. Indeed, the pro-cess becomes immediately much more complex when there are not only active andinactive states, but also interactions with other molecular species [4,5]. Comparingmultiple simulations, it can be easily spotted that spikes between different repeti-tions are not synchronized, and the time intervals between spikes of the same signalare not constant. These fluctuations, typical of random processes, can be difficultlyassociated to specific reactions when in presence of a heterogeneous population ofagents undergoing different interactions. Nonetheless, formation of patterns andinformation of the complexity of a non-ergodic system can be extracted.The former type of isoforms suggests studying the signals composed mainly ofpositive spikes emerging from zero in terms of stochastic analysis of large fluctua-tions and noise also when in presence of external perturbations. Nonetheless, thegeneralization of the analysis to other molecules requires us to cope with the nomi-nal values which are different from zero and change in time during the simulationsdue to variation of the equilibrium condition. For these reasons, we have developeda method to deal with any types of signals produced by this ABM. Signals’ transformation
The agents generated by the ABM contain time dependent information relative totheir position, velocity, belonging compartment, molecular activation states, bound scolani et al.
Page 5 of 31 state and an internal clock. If each agent has n variables, and there are S agents,then the state of the system can be represented in a n × S dimensional phase space.In order to reduce the dimensionality of the system’s phase space, in this study weneglected all the information relative to the position and velocity and we consid-ered only the molecular activation state, the bound state and, when biologicallyrelevant, the cellular compartment. To be precise, these quantities are categorical;therefore, we can label them with discrete indices. The signals analysed are thesum of molecules with the same labelling. Let us call each of these quantities statevariables which have values y i ( t ), where i is the index for the specific state variableand y is the value dependent on the time t . Given the stochastic origin of the simu-lations, the resulting time series are fluctuating signals, which can be decomposedin two components: the trend, Y i ( t ), and the fluctuations around the trend, f i ( t ).We introduce a method to process the signal in order to obtain a positive definedtime series of the noisy component by using the time average, a filter that removesthe random fluctuations faster then a time lapse ∆ T : y i ( t ) = b A ∆ T y i ( t ) ∆ T →∞ −−−−−→ Y i ( t ) , where the over bar represents the time averaged signal and the operator is themoving average operator: b A T = 1∆ T Z t +∆ T t ( · ) dt ′ . This approach allowed us to sequentially cut out fluctuations at different and spe-cific time scales and to analyse the results also for a system in an out of equilibriumcondition perturbed by external mechanical signals sensed by the integrins on themembrane of the osteoblast.The periodic step-wise external perturbation changes the equilibrium state ac-cordingly to the frequency, and, for fast variations compared to the response of themolecules, the system remains out of equilibrium for the whole simulation while themolecules try to reach the equilibrium condition. The average of the signal, beingnon constant in time, is one of the reasons why the system is not ergodic. There-fore, we first remove the trend from each signals in order to obtain the fluctuationsaround the trend.The choice of using a time average operator instead of the standard ensembleaverage has been dictated by biological motivations. Mathematically, an ensembleaverage is intended as an average among different regions which are statisticallyindependent, while, biologically speaking, it could be reinterpreted as an averageof non communicating parts of the system. Therefore, there are no real means forthe molecules with a finite interaction range to perceive a global average signal.On the contrary, a time average represents the slow component of the dynamicscumulated into the disposition of the molecules which can affect the relaxation oflarge deviations. Therefore, for a biological complex system, a time average operatorseams physically more meaningful for comparing fluctuations against the trend ofthe system. scolani et al.
Page 6 of 31
Subordination theory and detection of events
In this work, we propose a numerical method to disentangle random processes calledsubordinated processes [32] into their respective parent process and directing pro-cess in order to analyse the latter in terms of patterns of recurrence of events [33].Let us show how to generate a subordinated process [34–39]. First, we generatea leading process, which can be deterministic or derived from a random distribu-tion X so that a trajectory x has values x ( t ⋆i ) at discrete positive times t ⋆i with t ⋆i − t ⋆i − = ∆ t ⋆ > ∀ i ∈ N . If we rescale the time such that ∆ t ⋆ = 1 unit oftime, then t ⋆i = i and the parent process is described by x ( t ⋆i ) = x ( i ) where thesteps x ( t ⋆i ) − x ( t ⋆i − ) are tossed from X . Second, we generate in a similar way adirecting process T from a random distribution with positive support such that thetossed τ i defines a monotonically increasing function t ( t ⋆i ) = P ij =0 τ j . The subor-dinated function is defined as the trajectory x ( t ). Generally t ⋆ is called the naturaltime and t is called the physical time, which is the macroscopic and experimentallyobservable time [40, 41].In many physical and biological systems the variation of a signal is delayed bythe presence of internal structures and processes that are not easily accessible toexperimental observations or are neglected [42]; therefore, in these cases, we ex-perimentally deal with signals showing extra layers of complexity. Retrieving theprocesses X and T does not represent an easy task especially for real or in silicobiological models. One of the limitations is the quantification of data available forstatistical analysis. Another limit is that the sampling time, generally, does notcoincide with the times t ⋆i ; Third, experimentally it may be not evident what X represents and if any perturbation in the signal should be considered relevant orjust noise. For these reasons, the disentangling of the subordinated process is notunique, and we preferred to address the problem numerically. Our approach requiresto identify critical events in the system under analysis.In our ABM the number of the molecules in specific states fluctuates due torandom microscopical interactions between the molecules. These fluctuations aredriven at meso/microscopical scale by interactions among specific single particles,but experimentally they are elusive and present many technical problems in termsof instrumental resolution or capability of tracking large number of single particles.Rapid variations in the number of molecules in specific states affect the dynamic ofthe model, and their generation-relaxation towards a changing equilibrium conditionis characteristic of the model itself. Hence, it is reasonable for this ABM to definecritical events as the large and abrupt fluctuations in each molecular specie. Inorder to detect the events for a molecule i , first, we filtered the time series, y i ,with a moving average over a time lapse τ on each state variable , such that y ′ i ( t ) = b A τ y i ( t ), then we applied a second time the operator b A τ to the variation ofeach molecular specie i from its estimated mean ∆ y ′ i ( t ) = y ′ i ( t ) − y i ( t ) the result ofwhich is y ′′ i ( t ) = b A τ ∆ y ′ i ( t ), where the intervals of time are constrained by τ > τ .The first moving average has been used to de-trend the time series from the signal(intended as the component with slower dynamics). The subsequent second movingaverage has been applied to determine the amplitude σ i ( t ) of the second momentof the non-local noise fluctuations ∆ y ′ i ( t ) central with respect to y ′′ i ( t ). The termnon-local is intended as non-local in time, and it includes all the times belongingto [ t − τ / , t + τ /
2] for each time t . scolani et al. Page 7 of 31
The time series ∆ y ′ i ( t ) fluctuate around zero. The Hilbert transform, defined asthe integral operator b H = P V Z ∞−∞ ( · ) dτπ ( t − τ ) , has been used on the de-trended signals to derive the symmetric envelope of the de-trended signal, y env ( t ) = (cid:12)(cid:12)(cid:12) b H y ′ i ( t ) (cid:12)(cid:12)(cid:12) , from which it has been subsequently detracted.The result, y ′′′ ( t ) = y ′ ( t ) − y env ( t ), is a positive defined time series with enhancedprominent peaks and dumped valleys close to the zero axis, fig. 2. For each com-ponent y ′′′ i of the signal, the peaks are detected as events E , if the signal minusone standard deviation of the fluctuations drops after each peak, p n ( t ), below thecorresponding preceding peak value: E = { p n ( t ) , ∃ t | y ′′′ i ( t ) − σ i ( t ) < y ′′′ i ( t ) } . The component y ′′′ i has regions of small values; the beginning and ending of thoseregions are detected as events if the subtraction of one standard deviation of thefluctuations to the signal drops below or increases above the zero. The WTDsbetween peaks or the time extent of peaks are considered. To avoid artefacts relatedto fixed bin size and to be able to compare distributions obtained with differentnumerical parameters (in general having different binning sizes), we use a KernelDensity Estimator (KDE) defined as KDE( t ) = w R t + wt K ( t − τ )( · ) dτ such that R ∞−∞ K ( t ) dt = 1, and w is the support of the operator. Signal post-processing
The directing process, defined above as one component of the subordinated signaland here numerically obtained by detection of critical events, provides informationon the recurrences of large fluctuations of the number of molecules rapidly reachingor departing from a given node of the network, fig. 2. The waiting time distribution(WTD) of the directing process depends on the chosen parameters; therefore, it isrepresentative of the specific dynamics of the model.In our model, all the reaction occurrence times have a maximum time limit, butnot all interactions at the level of single agents are ruled by a single time distribution.Indeed, binding of molecules is driven by Brownian motions, while dissociationand relaxation of molecules are uniformly distributed, or have a constant value.Furthermore, the external perturbation maintains the system in a out of equilibriumcondition.Consequently, the directing processes of many molecules largely depart from aclassical Poisson process, where the critical events obtained by numerical detectionwould be gamma distributed, and characterize by multimodality.In our approach the integral distribution is first derived by weighing each inter-time τ between two consecutive events, y ′′′ ( t ) and y ′′′ ( t + τ ), with the amplitude ofthe fluctuation y ′′′ ( t ). The area τ × y ′′′ ( t ) is a measure of the minimum cumulatedtime that molecules in a given state i , which have been involved in the same criticalevent, are going to spend in any state different from i . We can see it as the cumulated scolani et al. Page 8 of 31 dispersion time of molecules sharing the same kinematics by departing from thesame node of the network.The non normalized distribution gives the amount of the total number of moleculeswith the same averaged dispersion time found inside the entire duration of the sam-pled signal. We would like to stress that all the signals have been de-trended andpositively defined. The WTD, ψ ( τ ) is the probability density obtained by normaliz-ing the integral distribution to 1. In order to verify the repeatability of the results,for each set of parameters, we repeated 10 independent simulations. For each sim-ulation, we computed the distributions, and we used the 10 independent resultsto derive the corresponding punctual confidence interval, (notice that the top andbottom of confidence intervals around the WTDs are not probability functions).Given the large amount of data simulated and results, we restricted the analysesprevalently on integrins, the RAF-MEK-ERK module of the molecular pathwayand RUNX2. Our decision is also related to the importance of these proteins andtheir interactions in transduction of the mechanical forces from the membrane tothe nucleus of the osteoblast as confirmed in the literature. Simulations and sensitivity analysis
The generic ABM platform FLAME (Flexible Large-scale Agent Modelling Envi-ronment) [43–45], capable of generating C based parallelisable code executable ondifferent parallel hardware architectures, was used to create and simulate the model.FLAME implements a discrete time fixed sweep update scheme meaning, in the sim-ulation, the time advances at each iteration of a fixed amount defined as “time unit”in which all agents are updated only once in an internal defined order [31, 46]. Oursimulations start with 8890 agents evolving with a unit time step corresponding to1 second for a period of 9 · s (approximatively 24 h ). The initial condition forthe position and velocity of the agents has been chosen at random from a uniformdistribution. The mechanical load used to perturb the system is mimicked by aperiodic square wave function defined by its magnitude M and oscillation period P .The proportions of agents belonging to specific molecular species, at time t = 0,are constant among all the simulations performed, and their respective values areshown in tab. 1. The values of the parameters for the baseline simulations are shownin tab. 2. The simulations were run at the SHeffield Advanced Research Computercluster (SHARC) on machines with 2 x Intel Xeon E5-2630 CPUs and 64GB ofRAM. Each simulation required 8GB of RAM, 96 hours of execution time andapproximatively 500GB of hard disk space.The effect on the fluctuations of the molecular populations were investigatedby changing the magnitude M and oscillation period P of the external pertur-bation and the molecular complex dissociation and relaxation times T • → • for: T RAF act +MEK d → MEK act , T MEK act → MEK d , T MEK act +ERK d → ERK act , T MEK act +ERK d → MEK act , T ERK act → ERK d , T ERK act +RUNX2 d → ERK act (fig. 1(b)).All combinations for N M = 2 values of the perturbation’s magnitude, N P = 3values of the perturbation’s period, and N T = 5 values for each of the 6 dissocia-tion/relaxation times has been analysed, and the spanned values of the parametersare reported in tab. 3. In order to account for the repeatability of the ABMs, each scolani et al. Page 9 of 31 simulation has been repeated 10 times for a total of 1800 simulations. Analyses andother results derived by the simulations are stored in FIGSHARE/Google Drive(link).
Phenomenological model and constraints
In our biological molecular model, the number of simulated molecules inside anosteoblast is much smaller than in the real biological system. The reasons for sucha decision is that running ABMs requires large resources in terms of memory andcomputational costs. This represents an upper limit in the number of agents whichcan be simulated. Another limit is given by spurious behaviours introduced by finitesize and finite volume effects [47, 48] representing an inferior boundary conditionfor the number of agents in the model. In our case, we have an average increasingnumber of agents during the simulation, which limited us in the initial maximumnumber of molecules. Nevertheless, the model does not explicitly include any finitevolume exclusion among molecules. Furthermore, the proportions of some moleculessuch as integrins in a specific state on the cell membrane are strongly driven by theexternal perturbation independently from the total amount of integrins simulated;this phenomenological aspect overrides many finite size effects. Consequently, weare in power to state that we can rescale the concentrations of the molecules andmaintain the same averaged dynamic of the system by rescaling the interactionrange R inter of the molecules accordingly [48, 49]. Results and discussion
Waiting time distributions and spectrograms
Dissociation times
During interaction between MEK and RAF, MEK is phosphorylated and conse-quently activated; after dissociation of the compound, activated MEK propagatesthe mechanotransduction signal. In fig. 3, we found the fluctuations of the com-plex MEK-RAF presents two distinct dynamics depending on the dissociation time.For dissociation time around T RAF act +MEK d → MEK act = 10 s , the WTD is a bi-modal distribution with modes around the recurrence times τ = 4 s and τ = 9 s ,fig. 3(c). The distribution ψ ( τ ) is approximately 15% larger than in τ . When T RAF act +MEK d → MEK act is increased from 90 s to 1320 s , the dynamics of the fluctu-ations changes completely; the peak at τ becomes negligible and the distributionhas only one relevant maximum at τ = 9 s . At T RAF act +MEK d → MEK act = 1320 s the integrated distributions and the WTD show that a large portion of the in-tertimes’ fluctuations are close to the mode and the occurrence of critical eventsis more regular. For T RAF act +MEK d → MEK act = 10 s , we also see the number ofmolecules with same dispersion time is below 200 for any delays and any periodof the mechanical stimulation, while it is larger than 400 at τ for dissociationtimes T RAF act +MEK d → MEK act > s .From the analyses of the occurrence of critical events, all mRNA sequencesbound to ribosomes show similar dynamics. As illustrated in fig. 4 and fig. 5, for T MEK act +ERK d → ERK act = 1320 s , the WTDs of mRNAs bound to ribosomes are tri-modal and the corresponding peaks are at τ = 4 s , τ = 10 s and τ = 18 s , respec-tively. The shapes of the WTDs remain constant at all the epochs simulated with scolani et al. Page 10 of 31 the exception of the mode at larger recurrence time τ which tends to fluctuate dueto stochastic noise. These properties hold true for P equal to 1000 s and 5000 s (seefig. 4(a) and fig. 4(b)), but for constant mechanical loads, the WTDs loose their tri-modality and become definitely bimodal after the age t = 5000 s as shown in fig. 5(a)and fig. 5(b). Similarly, in fig. 5(c) the WTDs for T MEK act +ERK d → MEK act = 1320 s and P = 200000 s have three modes at τ = 4 s , τ = 10 s and τ = 18 s . After thesystem reaches the age of t = 5000 s , the mode at τ disappears while the othertwo modes remain at equal τ , fig. 5(d). In all other cases, the WTDs of mRNAsbound to ribosomes are bimodal. Systems under constant mechanical loads withlarge activation times of MEK or ERK have a dynamics characterized by bursts oftranslation occurring at intertimes smaller than 10 seconds. On the contrary, forperiodic mechanical loads, there is a larger probability that abrupt variations inmRNA translation are separated by intertimes twice as long.Free mRNAs have completely different WTDs compared to their bounded coun-terparts. Generally, all the free mRNA WTDs have modes at intertimes equal to τ = 5 s , τ = 41 s and τ = 55 s which are the same at all epochs for any ofthe dissociation/relaxation times considered and all values of P , M and dissocia-tion/relaxation times T • → • simulated, fig. 6. The intertime range between τ and τ is characterized by multiple peaks changing with the age of the system, fig. 6(a).Except for T MEK act +ERK d → ERK act < s , in the WTDs of free mRNAs there is astable peak around τ = 10 s and an oscillating tail which slowly decays till the nextstable mode at τ . Instead, for values of T MEK act +ERK d → ERK act smaller than 300 s ,there is no stable peak at τ = 10 s and decaying tail becomes a long fluctuatingplateau, see fig. 6(b). Consequently, faster activation of ERK induces a more focusedtranscription of mRNAs, while at larger values of T MEK act +ERK d → ERK act , there is awide set of intertimes between τ and τ with high probability which causes non-rhythmic rapid variations of mRNA production. It is important to emphasis thatconstancy of behaviours among all types of mRNA, which, in the ABM, are iden-tically interacting molecules and thus are dynamically indistinguishable, shows theconsistency of the fluctuation analyses and the robustness of the method.Reducing the values of T • → • causes an increase of molecular interactions and acorresponding effect on the number of fluctuations in the signals. Such phenomenoncan be seen in non normalized histogram of bounded molecular complexes likeRAF+MEK, fig. 7(a) and fig. 8(a). Instead, free molecules show a direct relationbetween the size of the fluctuations and the dissociation/relaxation value, fig. 7(b)-(c) and fig. 8(b)-(c). Nevertheless, we have observed that the tendency of higheramount of events for smaller value of T • → • does not always hold, see fig. 9(a)-(c). In such cases, after the initial increase of the number of events, there is sat-uration. Decreasing T • → • further, there is a variation in the dynamics, and thequantity of critical events becomes directly proportional to the value of the relax-ation/dissociation time. Mechanical perturbations
The functional form of the external perturbation is a periodic square wave whichis completely characterized by the amplitude M and the period P , while the phaseis set to zero, and the minimal magnitude is equal to 100 µP a , see tab. 2. In the scolani et al. Page 11 of 31 range of values of P investigated, two were much shorter than the total duration ofeach epoch, and in one case the period was set such that the stimulation appliedwas constant for the whole simulation. The external force has the role of activatingthe mechanotransduction cascade and the dynamic of the system which, otherwise,would remain trapped in its initial condition.For all the analysed epochs and the simulated periods, the results indicated thatthe perturbation functions as a switch. When the perturbation goes to its minimum,integrin activation approaches zero, and the other molecules in the pathway followthe same dynamics. For all the molecules which are downstream in comparison tothe integrins, the larger their network distance is, the larger their delays are infollowing the dynamics of the external perturbation (fig. 1(a)).The passage of the mechanical loads from high to low and vice versa is of theorder of 1 time unit (1s), and the duration P of each regime is at least 1000 s long(see tab. 2 and tab. 3). Indeed, all the WTDs analysed have a support shorter thanall P considered, and in fig. 3 - 5, 7 - 11, the recurrence of events does not require atime larger than 100 s .The paucity and sparsity of the events characterized by the high/low transitionsof the external perturbation do have a minor effect on the intertimes τ betweenmolecular events even in cases where the dynamics of large quantities of moleculessynchronize with the external signal. On the contrary, during the time period P ,molecules have the chance to form compounds or dissociate multiple times, andthese events, induced by the faster molecular dynamics, account for most of thestatistics shown in the WTDs.The increase of the magnitude M of the perturbation resulted in an increase of thenumber of molecules activated which can be clearly observed from the variations inthe non normalized histograms, but this does not always correspond to a significantchange in the distributions of the intertimes, see third row in fig. 10.For active integrins, increasing P of a factor 40 resulted in reduction of the vari-ability among independent repetitions of the simulations. In the first and secondrows of fig. 10, the WTDs with P = 5000 s present a tail with a large confidenceinterval, and at large τ , the shapes largely fluctuate among different epochs. In-stead, for P = 200000 s , the waiting times produce clear bimodal distributions withreduced confidence interval at large τ resulting in no tail and a stable shape duringdifferent ages. The reduction of noise on active integrin intertimes is strongly sen-sitive to the length of the period P , but not on the magnitude of the perturbation,see third row in fig. 10. We noticed the sensitivity on P of active integrin intertimesdoes not depend on the value of T RAF act +MEK d → MEK act , T MEK act +ERK d → MEK act and T MEK act → MEK d . Instead, when the ERK cycle is slowed down, corresponding to T MEK act +ERK d → ERK act , T ERK act → ERK d and T ERK act +RUNX2 d → ERK act larger than thebaseline values, the WTDs’ noisy tails disappear, and the distributions are thesame for any value of the perturbation period.Increasing the period P implies larger lapse of time in which the mechanotrans-duction is switched on. Consequently, a larger number of events is expected. If thedynamics of the system does not change sensibly during the oscillations of the per-turbation, then the area under the histograms increases with the increase of P ,while the WTDs do not change. scolani et al. Page 12 of 31
Ageing
The integral distribution and the WTD are quantities computed over a time period,named epoch, which are much larger than the support of the represented distribu-tions and smaller than the total length of the simulation. Because the system isin a non equilibrium condition, we cannot state that the WTD does not dependon the epoch; therefore, the WTD can be rewritten as ψ t ( τ ) where t is the initialtime of the epoch of fixed extension used for the computation. Larger epochs reducethe lack of statistics and noise, while smaller epochs allow us to compute ψ t ( τ ) bychoosing the age t from a larger range of available values. We tested different epochsizes, and we compromised for a period of 10000 seconds.Ageing is characteristic of all the interacting agents, even though it is not so easy toestimate the effects of the passage of time in the histograms and WTDs, due to largeconfidence intervals and noise therein. Indeed, for some parameters and molecules,both the histograms and WTDs show alternate crescendo/decrescendo patternsand irregular cyclical variations. Even if the system is maintained in an out ofequilibrium condition, the dynamics of the system responds and adapts accordinglyto a mechanical perturbation.Ageing is a direct consequence of the initial conditions; nevertheless, given thatthe system is perennially far from the equilibrium, ageing is prevalently due toa variation of regimes in the dynamics of the system. Therefore, the state of thesystem at t = 0 should not be seen as a static initial condition, but as if the systemreached such a state due to a specific long standing dynamics. Indeed, we havesimulated a case where the cell had been stimulated with two consecutive trains ofoscillating perturbations with different frequencies. The resulting histograms showthere was no ageing at the end of the first perturbation and that ageing appearedduring the second regime.Among all molecules analysed, transcribed factors are characterized by ageingwith a most evident and clear progression. From the integral distributions in fig. 11,we observe that, with the increasing of age t there is a reduction on the numberof abrupt fluctuations. A similar ageing effect appears in all the mRNAs bound toribosomes. On the other hand, in many of the cases analysed, the correspondingWTD does not show sufficient discrepancies between different epochs. This meansall the dispersion times decrease of the same amount as time passes, and the ageingeffect is reduced to a slowing down of the system dynamics. As a direct consequenceof the fixed amount of total ribosomes simulated in the cytoplasm, ribosomes intheir free state also present an ageing effect, but in this case, the dispersion timesincrease 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 intertime value τ ;consequently, their WTD does not preserve the shape, but changes from a trimodaldistribution to a bimodal distribution as time passes, fig. 5. Conclusions
The dynamics of various molecules involved in the bone cell mechanotransductionhas been addressed in the framework of ABM. Preliminary analyses on the resultingsignals showed the external perturbation functions as a switch, but did not suggestany further synchronization on time scales smaller than the time period of the scolani et al.
Page 13 of 31 external mechanical loads. The amount of molecules in time were driven by a slowerdynamics frequently interrupted by asynchronous large fluctuations. Such behaviouris a common characteristic of subordinated processes. The lack of synchronizationhas suggested retrieving the directing processes of the signals and analyse theirdistributions in function of different parameters.The method of disentangling the principal process and the directing process is notunique, and it depends on the modality used to detect the events. The approachdeveloped in this work permits dealing with systems in out of equilibrium conditionsresulting from a particular initial state or due to the external perturbations, and itcould be applied to different types of time series. This method has been used to studythe dynamics of osteoblasts’ mechanotransduction which connects the extra cellularproperties with the intracellular signals and includes a large set of diffusive chemicalreactions. By identifying the occurrence of an event with multiple reactions at themolecular scale happening at the same time, the proposed numerical techniqueallowed us to derive the directing process for each molecular class. Differently frommore traditional approaches, where noise has been adopted as a measure to quantifythe error around the expected values of time dependent signals, here, the abruptfluctuations departing from the trend have been explored and analysed in orderto extract hidden information characteristic to each class of molecule and to theprocess dynamics.The WTD of many classes of molecules presented multimodality which highlightthe presence of preferred intertimes in the occurrence of events. Changes in therelaxation or dissociation times between molecules in the RAF-ERK-MEK moduleshowed variations in WTDs of molecules belonging both to the classes near oneof the perturbed edges and to other modules of the mechanotransduction network.For example, under a 1000 s time period perturbation, the WTD of RAF boundto MEK was characterized by a bimodal distribution for low values of active MEKdissociation times, while it became unimodal at larger value of the same parameter.Many of the mRNAs’ WTDs were affected by ageing, a variation of the shapefrom a bimodal to trimodal distribution, showing the emergence of a new preferredintertime of occurrence of event. Changes in the WTDs’ shape are proved to besignificant by local error estimation of the distribution.Variation in the modes implies variations in the dispersion times between chainedevents of the molecular network. Hence, the WTD of a specific molecule is a dy-namics’ marker which can be used as a signature for pathological conditions helpingto design drugs capable of improving target selectivity, maximizing the interactionwith targets and antagonizing chained reactions. Our method of analysis is a suit-able tool to study in vitro experiments tailored to investigate the effects of drugswith different properties like mass, diffusion coefficients or hydrosolubility. Thanksto FLAME’s parallelization and the genericity of the developed method, aspectslike the Parathyroid hormone or inflammatory factors involved in bone remod-elling and disease could be easily explored. Furthermore, the analyses done couldbe easily extended to identify the principal processes of the molecules which de-scribe the dynamics of the molecular baseline occurring at discrete time units. Suchdynamics identify discrete variations which can be linked to overall survival andworsening/improved survival analysis of diseases. scolani et al. Page 14 of 31
Competing interests
The authors declare that they have no competing interests.
Acknowledgements
This study was funded by the Engineering and Physical Sciences Research Council (EPSRC, Frontier MultisimGrant: EP/K03877X/1)
Author details Department of Oncology & Metabolism, University of Sheffield, Sheffield, UK. Insigneo Intitute for in silico medicine, University of Sheffield, Sheffield, UK. Department of Mechanical Engineering, University of Sheffield,Sheffield, UK.
References
1. Yavropoulou, M.P., Yovos, J.G.: The molecular basis of bone mechanotransduction. J. Musculoskelet. NeuronalInteract. (3), 221 (2016)2. Chen, J.C., Castillo, A.B., Jacobs, C.R.: Chapter 20 - Cellular and Molecular Mechanotransduction in Bone.Osteoporosis (Fourth Edition), 453–475 (2013). doi:10.1016/B978-0-12-415853-5.00020-03. Yap, A.S., Duszyc, K., Viasnoff, V.: Mechanosensing and Mechanotransduction at Cell–Cell Junctions. ColdSpring Harbor Perspect. Biol. (8), 028761 (2018). doi:10.1101/cshperspect.a0287614. Friedl, P., Mayor, R.: Tuning Collective Cell Migration by Cell–Cell Junction Regulation. Cold Spring HarborPerspect. Biol. (4), 029199 (2017). doi:10.1101/cshperspect.a0291995. Harris, A.R., Jreij, P., Fletcher, D.A.: Mechanotransduction by the Actin Cytoskeleton: Converting MechanicalStimuli into Biochemical Signals. Annu. Rev. Biophys. (1), 617–631 (2016).doi:10.1146/annurev-biophys-070816-0335476. Alonso, J.L., Goldmann, W.H.: Cellular mechanotransduction. Aims Press (2016).doi:10.3934/biophy.2016.1.507. Martino, F., Perestrelo, A.R., Vinarsk´y, V., Pagliari, S., Forte, G.: Cellular Mechanotransduction: From Tensionto Function. Front. Physiol. (2018). doi:10.3389/fphys.2018.008248. Humphrey, J.D., Dufresne, E.R., Schwartz, M.A.: Mechanotransduction and extracellular matrix homeostasis.Nat. Rev. Mol. Cell Biol. (12), 802 (2014). doi:10.1038/nrm38969. Vincent, T.L.: Targeting mechanotransduction pathways in osteoarthritis: a focus on the pericellular matrix.Curr. Opin. Pharmacol. (3), 449–454 (2013). doi:10.1016/j.coph.2013.01.01010. Papachroni, K.K., Karatzas, D.N., Papavassiliou, K.A., Basdra, E.K., Papavassiliou, A.G.:Mechanotransduction in osteoblast regulation and bone disease. Trends Mol. Med. (5), 208–216 (2009).doi:10.1016/j.molmed.2009.03.00111. Bonewald, L.F.: The Amazing Osteocyte. J. Bone Miner. Res. (2), 229 (2011). doi:10.1002/jbmr.32012. Oh, S.-H., Kim, J.-W., Kim, Y., Lee, M.N., Kook, M.-S., Choi, E.Y., Im, S.-Y., Koh, J.-T.: The extracellularmatrix protein Edil3 stimulates osteoblast differentiation through the integrin α β (11), 0188749 (2017). doi:10.1371/journal.pone.018874913. Charras, G., Yap, A.S.: Tensile Forces and Mechanotransduction at Cell–Cell Junctions. Curr. Biol. (8),445–457 (2018). doi:10.1016/j.cub.2018.02.00314. Kim, C., Ye, F., Ginsberg, M.H.: Regulation of Integrin Activation. Annu. Rev. Cell Dev. Biol. (1), 321–345(2011). doi:10.1146/annurev-cellbio-100109-10410415. Zhu, J., Luo, B.-H., Xiao, T., Zhang, C., Nishida, N., Springer, T.A.: Structure of a Complete IntegrinEctodomain in a Physiologic Resting State and Activation and Deactivation by Applied Forces. Molecular Cell (6), 849–861 (2008). doi:10.1016/j.molcel.2008.11.01816. Buo, A.M., Stains, J.P.: Gap Junctional Regulation of Signal Transduction in Bone Cells. FEBS Lett. (8),1315 (2014). doi:10.1016/j.febslet.2014.01.02517. Thiel, A., Reumann, M.K., Boskey, A., Wischmann, J., von Eisenhart-Rothe, R., Mayer-Kuckuk, P.: Osteoblastmigration in vertebrate bone. Biol. Rev. (1), 350–363 (2018). doi:10.1111/brv.1234518. Sun, Z., Guo, S.S., F¨assler, R.: Integrin-mediated mechanotransduction. J Cell Biol (4), 445–456 (2016).doi:10.1083/jcb.20160903719. Wang, N.: Review of cellular mechanotransduction. J. Phys. D: Appl. Phys. (23), 233002 (2017).doi:10.1088/1361-6463/aa6e1820. Muhamed, I., Chowdhury, F., Maruthamuthu, V.: Biophysical Tools to Study Cellular Mechanotransduction.Bioengineering (1), 12 (2017). doi:10.3390/bioengineering401001221. Shams, H., Soheilypour, M., Peyro, M., Moussavi-Baygi, R., Mofrad, M.R.K.: Looking “Under the Hood” ofCellular Mechanotransduction with Computational Tools: A Systems Biomechanics Approach across MultipleScales. ACS Biomater. Sci. Eng. (11), 2712–2726 (2017). doi:10.1021/acsbiomaterials.7b0011722. D’Antonio, G., Macklin, P., Preziosi, L.: An agent-based model for elasto-plastic mechanical interactionsbetween cells, basement membrane and extracellular matrix. Math. Biosci. Eng. (1), 75–101 (2013).doi:10.3934/mbe.2013.10.7523. Piel, M., Thery, M.: Methods in cell biology micropatterning in cell biology part b volume 120 preface. Methodsin cell biology , (2014). doi:10.1016/B978-0-12-417136-7.09989-424. Jurchenko, C., Salaita, K.S.: Lighting up the force: Investigating mechanisms of mechanotransduction usingfluorescent tension probes. Mol. Cell. Biol. (15), 2570–2582 (2015). doi:10.1128/MCB.00195-1525. Kamkin, A., Kiseleva, I.: Mechanosensitivity and Mechanotransduction. Springer, (2011).doi:10.1007/978-90-481-9881-826. Peng, T., Liu, L., MacLean, A.L., Wong, C.W., Zhao, W., Nie, Q.: A mathematical model ofmechanotransduction reveals how mechanical memory regulates mesenchymal stem cell fate decisions. BMCSyst. Biol. (2017). doi:10.1186/s12918-017-0429-x scolani et al. Page 15 of 31
27. Dong, X., Foteinou, P.T., Calvano, S.E., Lowry, S.F., Androulakis, I.P.: Agent-Based Modeling ofEndotoxin-Induced Acute Inflammatory Response in Human Blood Leukocytes. PLoS One (2), 9249 (2010).doi:10.1371/journal.pone.000924928. Ausk, B.J., Gross, T.S., Srinivasan, S.: An agent based model for real-time signaling induced in osteocyticnetworks by mechanical stimuli. J. Biomech. (14), 2638–2646 (2006). doi:10.1016/j.jbiomech.2005.08.02329. Zhang, L., Wang, Z., Sagotsky, J.A., Deisboeck, T.S.: Multiscale agent-based cancer modeling. J. Math. Biol. (4-5), 545–559 (2009). doi:10.1007/s00285-008-0211-130. Bonchev, D., Thomas, S., Apte, A., Kier, L.B.: Cellular automata modelling of biomolecular networksdynamics. SAR QSAR Environ. Res. (1-2), 77–102 (2010). doi:10.1080/1062936090356858031. FLAME-HPC github repository. https://github.com/FLAME-HPC/xparser. Comments and discussions inhttps://github.com/svdhoog/xparser32. Sokolov, I.M., Klafter, J.: From diffusion to anomalous diffusion: A century after einstein’s brownian motion.Chaos: An Interdisciplinary Journal of Nonlinear Science (2), 026103 (2005). doi:10.1063/1.1860472.https://doi.org/10.1063/1.186047233. Gorenflo, R., Mainardi, F.: Subordination pathways to fractional diffusion. Eur. Phys. J. Spec. Top. (1),119–132 (2011). doi:10.1140/epjst/e2011-01386-234. Montroll, E.W., Weiss, G.H.: Random walks on lattices. ii. Journal of Mathematical Physics , 1174 (2018). doi:10.3389/fphys.2018.0117439. Ascolani, G., Bologna, M., Grigolini, P.: Subordination to periodic processes and synchronization. Physica A:Statistical Mechanics and its Applications , 7366 (2015). Article43. FLAME. flame.ac.uk44. Bajo, J., Corchado, J.M., Mart´ınez, E.M.N., Icedo, E.O., Mathieu, P., Hoffa-Dabrowska, P., Val, E., Giroux, S.,Castro, A.J.M., S´anchez-Pi, N., et al. : Highlights of Practical Applications of Agents, Multi-Agent Systems,and Complexity: The PAAMS Collection: International Workshops of PAAMS 2018, Toledo, Spain, June 20–22,2018, Proceedings. Communications in Computer and Information Science. Springer, (2018).https://books.google.co.uk/books?id=KutgDwAAQBAJ45. Richmond, P., Walker, D., Coakley, S., Romano, D.: High performance cellular level agent-based simulationwith flame for the gpu. Briefings in Bioinformatics (3), 334–347 (2010). doi:10.1093/bib/bbp07346. Fat`es, N.: A guided tour of asynchronous cellular automata. In: Kari, J., Kutrib, M., Malcher, A. (eds.) CellularAutomata and Discrete Complex Systems, pp. 15–30. Springer, Berlin, Heidelberg (2013)47. Ascolani, G., Badoual, M., Deroulers, C.: Exclusion processes: Short-range correlations induced by adhesion andcontact interactions. Phys. Rev. E , 012702 (2013). doi:10.1103/PhysRevE.87.01270248. Rhodes, D.M., Holcombe, M., Qwarnstrom, E.E.: Reducing complexity in an agent based reactionmodel—benefits and limitations of simplifications in relation to run time and system level output. Biosystems , 21–27 (2016). doi:10.1016/j.biosystems.2016.06.00249. Pogson, M., Smallwood, R., Qwarnstrom, E., Holcombe, M.: Formal agent-based modelling of intracellularchemical interactions. Biosystems (1), 37–45 (2006). doi:10.1016/j.biosystems.2006.02.004. Dedicated to theMemory of Ray Paton Figures scolani et al.
Page 16 of 31(a)(b)
RAF act + MEK d MEK act
MEK act + ERK d MEK d ERK d ERK act
ERK act + RUNX2 d T RAF act +MEK d → MEK act T MEK act → MEK d T MEK act +ERK d → ERK act T MEK act +ERK d → MEK act T ERK act → ERK d T ERK act +RUNX2 d → ERK act
Figure 1 Mechanotransudction molecular network.
Schematic representation of all the molecularspecies simulated in the ABM (white boxes). The arrows show the interactions. Multiple blackarrows entering in a white boxes describes mass action reactions, while red arrows coming out (orarrowheads) means dissociation of molecular complexes or energetic relaxation of the molecule.Orange lines separates osteblast compartments. In the coloured boxes the network modules areenclosed. (a)
Complete diagram with the external stimuli (cyan), the energetic MAPK mod.(pink), the transcription mod. (blu), the translation mod. (green) and the ECM proteins (red). (b)
Part of the diagram showing only the simulation dependent delays and the respectivemolecules. Interactions depicted in black have fixed parameters.scolani et al.
Page 17 of 31 ph ys i c a l t i m e - t [ s ] M E K bound E RK s ubo r d i n a t e d p r o cess t r end + en v e l ope o f m o l e c u l e s ( s ubo r d i na t ed p r o c e ss ) t r end na t u r a l t i m e - t * i * i )& size of fluctuations p r i n c i p a l p r o cess & d i r ec t i ng p r o cess s t e p s t r end ( p r i n c i pa l p r o c e ss ) de t r ended o f m o l e c u l e s i n t e r t i m e s be t w een ab r up t pea ks ( d i r e c t i ng p r o c e ss s t ep s ) intertimes - i =t i -t i-1 [s] V i s ua l app r o x i m . o f s t r e t c h i ng / s h r i n k i ngo f t he s i gna l due t o i n t e r t i m e s ' ab l a t i on Figure 2 Subordinated process.
On the top graph, the amount of molecules for the complexMEK+ERK in function of the time t is shown . This curve represents the subordinated signal x ( t ) . The trend of the signal, derived by time average, and the envelope, derived by Hilberttransform, are used to detect the critical events. In the bottom graph, the intertimes τ i betweentwo consecutive critical events are shown in function of the physical time t ⋆i = i . The curve t ( t ⋆i ) ,given by the sum of all the intertimes τ up to the i-th event, is a realization of the directingprocess. In the same graph, the curve x ( t ⋆i ) , given by the sum of all the steps ∆ x up to the i-thevent, is a realization of the principal process. For convenience, in the bottom graph, thefluctuations around the trend are shown. In the curve x ( t ⋆i ) , regions with high number of eventsresemble a stretched form of the subordinated signal while regions with few events resembleshrunk portions of the subordinated signal.scolani et al. Page 18 of 31(a)(b)(c)(d)
Figure 3 MEK bound RAF.
Independently from the frequencies of the tested perturbations, thevariation of T RAF act +MEK d → MEK act leaves the MAPK sub module of the network unchangedexception for the MEK interacting with RAF which is highly sensitive to the parameter. Indeed,the distribution of recurrence of events for active MEK bounded to RAF is bimodal in τ = { , } s . For small value of T RAF act +MEK d → MEK act ∼ s the distribution is mainlycentred around τ = 4 s and it represents the main mode, while for larger value of T RAF act +MEK d → MEK act the distribution’s main mode is at τ = 9 s . (a) T = 1000 s ; (b) T = 200000 s ; (c) T RAF act +MEK d → MEK act = 10 s ; and (d) T RAF act +MEK d → MEK act = 1320 s .scolani et al. Page 19 of 31(a)(b)
Figure 4 Ageing of ALP mRNA.
For T MEK act +ERK d → ERK act = 1320 s , the WTDs of the ALPmRNA have trimodal distributions with peaks at s , s and s . The three modes arepersistent at any age of the system. For T MEK act +ERK d → ERK act ≤ s , the WTDs’ slope at τ > s remains negative and the shape of the WTDs is bimodal. (a) and (b) M = 10000 µP a and P = 1000 s , while T ERK act → ERK d changes; (b) aged condition of the case on the above line.scolani et al. Page 20 of 31(a)(b)(c)(d)
Figure 5 Ageing of ALP mRNA.
For T MEK act +ERK d → MEK act = T MEK act +ERK d → ERK act = 1320 s , the WTDs of the ALP mRNA at t = 0 s have trimodaldistributions with peaks at s , s and s . After the system reaches age t = 5000 s , the WTDs’slope at τ > s remains negative. In the other cases, the peak at τ = 18 s is not present. (a) - (d) M = 10000 µP a and P = 200000 ; (a) T MEK act +ERK d → ERK act changes; (b) aged condition ofthe case on the above line; (c) T MEK act +ERK d → MEK act changes; (d) aged condition of the caseon the above line.scolani et al.
Page 21 of 31(a)(b)
Figure 6 BSP mRNA.
The BSP mRNA WTDs have 4 modes. The shape of the distributions aresimilar to all the unbounded mRNA in the cytoplasm. The modes are at τ equal to s , s , s and s . In the region between the modes at s and s , the distributions are noisy and theiraverage slope depend on T • → • . (a) P = 1000 s , M = 10000 µP a and T MEK act → MEK d change; (b) P = 200000 s , M = 10000 µP a and T MEK act +ERK d → ERK act change.scolani et al.
Page 22 of 31(a)(b)(c)
Figure 7 MEK+ERK, ERK and RUNX2.
Distributions for different values of T ERK act → ERK d under a periodic perturbation with P = 1000 s and M = 10000 µP a . T ERK act → ERK d produces larger effects on the next nearest neighbour of the network likeunphosphorylated RUNX2 and unphosphorylated ERK. Variations in the distributions and in thenon-normalized histogram of MEK bounded to ERK. No effects are visible on the distribution ofactive ERK dissociated from MEK. Amplitudes or number of fluctuations of ERK dissociatedMEK increase with the increase of T ERK act → ERK d , while they show an inverse relation for MEKbound to ERK decrease. (a) MEK bound to
ERK ; (b) active ERK dissociated from
MEK ; (c) unphosphorylated RUNX2 .scolani et al.
Page 23 of 31(a)(b)(c)
Figure 8 MEK+ERK, ERK and RUNX2.
Distributions for different values of T ERK act → ERK d under a periodic perturbation with P = 200000 s and M = 10000 µP a . Thedistributions are similar to those obtained with higher frequency stimulation (see fig. 7). Thehistograms show larger or more frequent fluctuations of unphosphorylated RUNX2 and MEKbound ERK. (a) MEK bound to
ERK ; (b) active ERK dissociated from
MEK ; (c) unphosphorylated RUNX2 .scolani et al.
Page 24 of 31(a)(b)(c)
Figure 9 ERK+MEK complex.
The distributions of intertimes for the ERK+MEK complexWTDs are bimodal. The modes are at τ equal to s and s . When T MEK act +ERK d → ERK act = 8 s , the mode at τ = 15 s shifts to s . The non normalizedhistograms on the left column show that the total number of events is at the minimum when T MEK act +ERK d → ERK act has the maximum value simulated equal to s . The number ofcritical events increases when T MEK act +ERK d → ERK act decreases and reaches a critical value afterwhich the dissociation time is directly proportional to the number of events. (a) - (c) M = 10000 µP a ; (a) P = 1000 s ; (b) P = 5000 s ; (c) P = 200000 s .scolani et al. Page 25 of 31(a)(b)(c)
Figure 10 Active integrins.
The distribution of recurrence of events for active integrins is bimodalin τ ∼ { , } s . The modes are independent from the periods and the magnitude of themechanical load. The tail of the distributions at large τ are sensitive to the perturbation period P .For P = 5000 s , the distributions show a noisy large tail. For P = 200000 s , the distributions decayrapidly at zero around tau ∼ s . At large τ , the distribution does not depend on the magnitude M . (a) M = 1000 µP a while P and T MEK act +ERK d → MEK act change; (b) M = 1000 µP a while P and T MEK act → MEK d change; (c) P = 200000 s while M and T MEK act → MEK d change.scolani et al. Page 26 of 31
Figure 11 BSP mRNA bound to ribosome.
Spectrogram of the critical events for BSP mRNAduring translation. On the x axis, the age t of the non normalized histogram also corresponding tothe beginning of an epoch; on the y axis, the intertimes τ ; on the z axis, the number of eventsregistered during each epoch of length s . The vertical black lines crossing the surface showthe 95% CI. The color defines the value of the WTD ψ at intertime τ and age t . The WTD isbimodal and the shape is constant at all t . Nonetheless, the total number of critical eventsdecreases with the age.scolani et al. Page 27 of 31
Tables
Table 1
Number of molecules at time t = 0 .Molecule Osteoblast cell I d I f act I f act + FAK act I act FAK act I f act + FAK act I f act + FAK act + RAS d FAK d RAS d I f act + FAK act + RAS d RAS act RAS act + RAF d RAF d RAS act + RAF d RAF act RAF act + MEK d MEK d RAF act + MEK d MEK act MEK act + ERK d ERK d MEK act + ERK d ERK act ERK act + RUNX2 d ERK act + OTHER RUNX2 d ERK act + RUNX2 d RUNX2 act RUNX2 act + DNA OSX mon RUNX2 act + OSX mon OSX mon + DNA RUNX2 act + OSX mon + DNA RUNX2 act + OSX mon + DNA + AP1 OSX mul RUNX2 act + OSX mul OSX mul + DNA RUNX2 act + OSX mul + DNA AP1 RUNX2 act + OSX mon + DNA + AP1 OPN
RNA Rib av + OPN RNA Rib n av + OPN
RNA OCN
RNA Rib av + OCN RNA Rib n av + OCN
RNA ALP
RNA Rib av + ALP RNA Rib n av + ALP
RNA BSP
RNA Rib av + BSP RNA Rib n av + BSP
RNA Rib n comp
Rib n comp + RNA Rib c Rib c + RNA MAT + Vesc et al.
Page 28 of 31
Table 2
Parameters’ names, symbols and values.Parameter Symbol Value Unit ofmeasureMechanical min. magnitude m µP a
Mechanical phase φ T I d → I f act I d dissoc. time from I + FAK comp. T I f act +FAK act → I d U [0 , s FAK act activation delay T FAK d → FAK act
FAK d dissoc. time from I + FAK comp. T I f act +FAK act → FAK d
30 s I act + FAK dissoc. time T I f act +FAK act +RAS d → I f act +FAK act RAS act dissoc. time from
I+ FAK+ RAS comp. T I f act +FAK act +RAS d → RAS act
RAS act dissoc. form
RAS + RAF comp. T RAS act +RAF d → RAS act
14 s
RAS act relaxation time T RAS act → RAS d
60 s
RAF act dissoc. time T RAS act +RAF d → RAF act
14 s
RAF act dissoc. time from
RAF + MEK T RAF act +MEK d → RAF act
10 sRelaxation of
RAF act T RAF act → RAF d
60 s
MEK act dissoc. time from
RAF + MEK comp. T RAF act +MEK d → MEK act
10 s
MEK act dissoc. time from
MEK + ERK comp. T MEK act +ERK d → MEK act
MEK relaxation time T MEK act → MEK d
88 s
ERK activation time T MEK act +ERK d → ERK act
ERK relaxation time T ERK act → ERK d
50 s
ERK d dissoc. time from RUNX2 + ERK comp. T ERK act +RUNX2 d → ERK act
10 s
ERK act dissoc. time from
ERK act + OTHER comp. T ERK act +OTHER → ERK act
RUNX2 act dissoc. time from
ERK act + RUNX2 d comp. T ERK act +RUNX2 d → RUNX2 act
10 s
RUNX2 act dissoc. time from
RUNX2 act + DNA comp. T RUNX2 act +DNA → RUNX2 act
10 s
RUNX2 relaxation time T RUNX2 act → RUNX2 d
60 s
OSX mon dissoc. time from
RUNX2 act + OSX mon comp. T RUNX2 act +OSX mon → OSX mon
20 s
OSX mon + DNA dissoc. time from
RUNX2 act T RUNX2 act +OSX mon +DNA → OSX mon +DNA
30 s
OSX mon + DNA dissoc. time from
AP1 T RUNX2 act +OSX mon +DNA+AP1 → OSX mon +DNA
OSX mul dissoc. time from
RUNX2 act + OSX mul comp. T RUNX2 act +OSX mul → OSX mul
40 s
OSX mul + DNA dissoc. time from
RUNX2 act T RUNX2 act +OSX mul +DNA → OSX mul +DNA
50 sribosome availability time from
Rib n av + OPN
RNA T Rib n av +OPN
RNA → Rib av +OPN RNA
20 s
OPN
RNA dissociation time from available ribosome T Rib av +OPN RNA → OPN
RNA
75 sribosome availability time from
Rib n av + OCN
RNA T Rib n av +OCN
RNA → Rib av +OCN RNA
20 s
OCN
RNA dissociation time from available ribosome T Rib av +OCN RNA → OCN
RNA
75 sribosome availability time from
Rib n av + ALP
RNA T Rib n av +ALP
RNA → Rib av +ALP RNA
20 s
ALP
RNA dissociation time from available ribosome T Rib av +ALP RNA → ALP
RNA
75 sribosome availability time from
Rib n av + BSP
RNA T Rib n av +BSP
RNA → Rib av +BSP RNA
20 s
BSP
RNA dissociation time from available ribosome T Rib av +BSP RNA → BSP
RNA
75 s
RNA dissociation time from
Rib n comp T Rib n comp +RNA → Rib n comp
100 s
RNA dissociation time from complete ribosome T Rib c +RNA → Rib c
100 sforce tag probability F tag P { , } = { . , . } I d interaction radius R I d
50 nm
FAK act interaction radius R FAK act
50 nm
RAS d interaction radius R RAS d
80 nm
RAF d interaction radius R RAF d
70 nm
MEK d interaction radius R MEK d
50 nm
ERK d interaction radius R ERK d
30 nm
RUNX2 d interaction radius R RUNX2 d
30 nm
OSX mon interaction radius R OSX mon
30 nm
Rib n comp interaction radius R Rib n comp
30 nmcell radius R cell R nucl
400 nmprotein’ s average velocity v p ∆ v p ∆ ϕ p = ∆ θ p π /10 radprotein translation delay T p
95 stranscription delay per mRNA T RNA
600 sIntegrins θ angle on cell membrane distribution I θ U [0 , π ] φ angle on cell membrane distribution I φ U [0 , π ] et al. Page 29 of 31
Table 3
Parameters ranges: Names, symbols, unit of measures and list of values simulated. Boldquantities represent the baseline values. Where no baseline is present, then all possible combinationshas been considered. Each set of parameters has been independently repeated 10 times.Parameter Symbol List of values Unit ofmeasureMechanical max. magnitude M { , } µP a Mechanical period P { , , } s MEK act dissoc. time from
RAF act + MEK d comp. T RAF act +MEK d → MEK act { , , , , } s MEK act dissoc. time from
MEK act + ERK d comp. T MEK act +ERK d → MEK act { , , , , } s MEK d relaxation time T MEK act → MEK d { , , , , } s ERK act activation time T MEK act +ERK d → ERK act { , , , , } s ERK d relaxation time T ERK act → ERK d { , , , , } s ERK act dissoc. time from
ERK act + RUNX2 d comp. T ERK act +RUNX2 d → ERK act { , , , , } sscolani et al. Page 30 of 31 I d I act I f act I f act + FAK act FAK act I f act + FAK act I f act + FAK act + RAS d FAK d RAS d I f act + FAK act + RAS d RAS act RAS act + RAF d RAF d RAS act + RAF d RAF act RAF act + MEK d MEK d RAF act + MEK d MEK act MEK act + ERK d ERK d MEK act + ERK d ERK act ERK act + RUNX2 d ERK act + OTHER RUNX2 d ERK act + RUNX2 d RUNX2 act RUNX2 act + DNA act actRD and R p-R AA R p-A p0 D and R p-D p0substrate p+ D and R p-A D p+R p- AD p0D p+D p-R p-A AA D p+R p- AD p0D p+ A p0R p- D p0 D p0A AD p+R p- DA A Figure 12
Reactions in the ECM and cytoplasm. A=association, D=dissociation, R=relaxation,p+=phosphorylation, p-=dephosphorylation, p0=no phosphorylation, +1=creation,act=mechanotransduction activation, rand=chosen at randomscolani et al.
Page 31 of 31
RUNX2 act + OSX mon
RUNX2 act + OSX mon + DNA
RUNX2 act + OSX mul
OSX mon OSX mon + DNA OSX mul OSX mul + DNA RUNX2 act + OSX mon + DNA + AP1
RUNX2 act + OSX mul + DNA
OPN
RNA Rib av + OPN RNA Rib n av + OPN
RNA OCN
RNA
Rib av + OCN RNA
Rib n av + OCN
RNA
ALP
RNA
Rib av + ALP RNA
Rib n av + ALP
RNA
BSP
RNA
Rib av + BSP RNA
Rib n av + BSP
RNA
Rib n comp
Rib n comp + RNA
Rib c Rib c + RNA MAT + Vesc
D D D DD D D DA A A Arand and +1DD A ADD A ADD A ADD A AD DA A +