Search

Additional file 16: Fig. S7. Pseudo-time evolution of variables during simulations. The sequence of variable updating over each time steps after introducing a perturbation from the healthy state was saved as a timeseries and plotted. It shows the discrete behavior of the simulation and that a steady state is reached way before the duration of the perturbation (1000 time steps) is reached in such a way that maintain the perturbation longer would not change the output of the simulation once the perturbation is released.

Additional file 12: Table S4. Targets and associated small molecules or growth factors for in vitro validation. The first row indicates the targets to be perturbed, [+] stands for activation while [-] stands for inhibition. The name (resp. cat number) of the small molecule or growth factor employed to achieve that effect is indicated in the row called ‘Molecule name’ (resp. ‘Cat n°’).

Additional file 3: Table S2. List of variables and mouse genes correspondence. All mathematical variables and the corresponding node in the network have names written in upper cases and do not reflect the official human or mouse nomenclatures. To relate those variables to actual genes more easily, we provide this table of correspondence. A related mouse gene name and NCBI ID is indicated for each variable. Nevertheless, this is not exhaustive since some variables represent a group of factors or a family of ligands rather than a single factor.

Additional file 7: Table S3. Validation of transcriptional data-inferred regulatory interactions. Transcriptional interactions were inferred from the merged OA dataset. Inference was run with three algorithms and only interactions that were present in the results of the three algorithms (GENIE3, ARACNE, TIGRESS) were considered as additions to the model. An interaction was considered present for one algorithm if it scored higher than a threshold defined as the difference between the mean and standard deviation of all scores (Additional file 6, Data S2 ). Inferred interactions were validated whenever possible by looking for binding sites of the source transcription factor in the enhancer region of the target gene with GeneCards. The complete list of GeneHancer IDs provided by GeneCards are reported in the last column. ‘NA’ indicate that no binding site was found. The exact gene names and gene IDs that were queried in GeneCards are reported in the first columns.

Additional file 10: Fig. S2. Effect of ALKs ratio in the influence of inflammation and TGFβ signaling on chondrocyte hypertrophy. Percentage of perturbations remaining in the healthy state (right) or transitioning towards the hypertrophic one (left) during inflammatory pathway activation with TGF-B treatment while changing the ratio between ALK1 and ALK5. The inflammatory and TGF-B profiles that were imposed are the same as in Fig. 4A except that the value imposed for ALK1 and ALK5 are varying between 0 and 1 with a 0.1 increment. ALK1=ALK5 on the red diagonals and ALK1>ALK5 in the upper right corner. Roughly, the rescue of the healthy state by TGF-B is lost when ALK1 is greater than ALK5. If ALK1 is high enough and that the difference between ALK1 and ALK5 is not greater than 10-20% then the protective effect of TGF-B is also disturbed for ALK1

Additional file 5: Fig. S1. Individual dataset clustering and heatmap. The microarray sub-datasets, used in the merged dataset, were also investigated individually. They were pre-processed (i.e. processing steps before merging and correcting for batch effect) in the same way as the merged dataset but the unsupervised clustering analysis was done on each sub-dataset separately in order to compare the biological or OA related information content before and after the data merging. The heatmaps show the expression profile of the same list of genes of interest than for the merged dataset (see Additional file 4, Data S1.). These expression datasets were submitted to unsupervised clustering with the Euclidean distance method and the Complete aggregation method in R thanks to the heatmap3 function from the Github repository https://github.com/obigriffith/biostar-tutorials/tree/master/Heatmaps . The headers indicate the GEO accession numbers of the 6 original datasets. Samples labeled as ‘WT’, for wild type, are in green, samples labeled as ‘OA’, for osteoarthritis, in red. The pre-labelling is the same as for the merged dataset (see Additional file 4, Data S1). When applicable, a grey scale indicate the time points (w stands for weeks, in the legend). For some datasets (e.g. GSE26475, GSE33656, GSE53857) the OA and the WT samples are well separated in different clusters, while for other the separation is not so clear. For instance, in GSE45793 the variance due the time point 6weeks is greater than the OA induced variance, while for weeks 1 and 2 OA and WT samples are well separated due to the OA condition.

Additional file 2. Supplementary computational method. Equations, mathematical framework and justification of deviation from the general rule in the equations.

Additional file 11: Fig. S3. Screenshot of the user-friendly interface for the virtual chondrocytes App. The standalone Matlab-based applications can be launched and used without Matlab license, provided that the compiler Matlab Runtime is installed ( https://nl.mathworks.com/products/compiler/matlab-runtime.html ). The virtual chondrocyte initial state can be set as healthy or hypertrophic, allowing the user to test any scenarios. All the 60 components may be perturbed alone or in any sort of combination by forcing the variables to take a value in the interval [0:1], with a step of 0.1. The most left column indicate the value of the variable in the selected initial state, for information. Obviously, applying a perturbation that is equal to the initial value of the variable will not affect the system. Once the setting are done, the user can apply the experimental condition by pushing the button ‘Test condition’ and the percentage of transitions towards each of the possible basal stable states (i.e. ‘None’, ‘Healthy’ and ‘Hypertrophic’) is computed. If the ‘Compute statistics’ box is ticked, then the experiment is repeated 3 times and the average and standard deviations are displayed (variation occurs due to the stochastic nature of the model). The results may be exported and saved in an excel file via the ‘Save’ button. The application can be installed with the executable file on the GitHub repository, [ https://github.com/Rapha-L/Insilico_chondro.git ]. No Matlab license is required, however, the operating system should be able to support the Matlab software. For Linux users, a Windows virtual machine may be used.

Additional file 5: Fig. S1. Individual dataset clustering and heatmap. The microarray sub-datasets, used in the merged dataset, were also investigated individually. They were pre-processed (i.e. processing steps before merging and correcting for batch effect) in the same way as the merged dataset but the unsupervised clustering analysis was done on each sub-dataset separately in order to compare the biological or OA related information content before and after the data merging. The heatmaps show the expression profile of the same list of genes of interest than for the merged dataset (see Additional file 4, Data S1.). These expression datasets were submitted to unsupervised clustering with the Euclidean distance method and the Complete aggregation method in R thanks to the heatmap3 function from the Github repository https://github.com/obigriffith/biostar-tutorials/tree/master/Heatmaps . The headers indicate the GEO accession numbers of the 6 original datasets. Samples labeled as ‘WT’, for wild type, are in green, samples labeled as ‘OA’, for osteoarthritis, in red. The pre-labelling is the same as for the merged dataset (see Additional file 4, Data S1). When applicable, a grey scale indicate the time points (w stands for weeks, in the legend). For some datasets (e.g. GSE26475, GSE33656, GSE53857) the OA and the WT samples are well separated in different clusters, while for other the separation is not so clear. For instance, in GSE45793 the variance due the time point 6weeks is greater than the OA induced variance, while for weeks 1 and 2 OA and WT samples are well separated due to the OA condition.

Additional file 1: Table S1. Definitions. Definitions of technical terms pertaining to biological signaling network modeling with a semi-quantitative additive formalism.

Additional file 14: Fig. S5. Illustration of the algorithm for the asynchronous updating of variables with a simplified (3 nodes) example network. The network represents interactions happening in one of the subnetworks (protein= fast reactions or genetic = slow reactions). Inhibitions are represented in red and activations are in black. The mathematical rules corresponding to the network are displayed. If the rules result in a value lower than 0 (resp. higher than 1), the value is brought back to 0 (resp. 1). For this example, 3 different initial states are inputted and each of the three variables is updated asynchronously. The order in which variables are updated is random. The system reaches a stable state when the next update gives the same state as in the previous time step and that all variables were screened in the random ordering list. That state is a pseudo-stable state if the rules were describing fast reactions, in that case, a new slow variable can be updated (see Fig. S4.). However it is a final stable state if the rules were describing slow interactions since it would mean that the system had first reached a pseudo-stable state at the fast level and would now be stable at the slow level too. The example illustrates that different initial states may reach the same final state but also that the same initial state (e.g. [1 0 1] ) can reach different final states, depending on the order in which variables are updated, thereby introducing stochasticity in the system.

Additional file 13: Fig. S4. Decision tree summarizing the variable updating scheme employed in algorithm to simulate the in silico chondrocyte. Each biological component is represented by the gene expression level (slow variable) and the protein activity potential (fast variable). Variables are updated based on the rules stored in the model’s equations. First, fast variables are updated in random order, when a pseudo-stable state is reached and that all fast variables have been updated, the next random chosen slow variable is updated. This goes on until a state that is stable both at the fast and slow level is reached. This is the final stable state. A state is considered stable if further variable updates do not bring further changes for any of the variables, with a predefined tolerance interval. The order in which variables are updated is random, thereby generating some stochasticity in the model. Within the fast (resp. slow) updating loops, variables are updated asynchronously (meaning the one after the others) according to the rules defined in the system of equations and in a random order. For some systems (i.e. set of equations) cyclic attractors may arise, meaning that the system never reaches a fixed stable state but oscillates between several states. This situation did not occur in the current study.

Additional file 15: Fig. S6. Sensitivity analysis: impact of the number of initialization on the canalization results during the Monte Carlo. The percentage of random initialization reaching each attractor is displayed with the Sox9 positive state (healthy) in orange, the Runx2 positive state in blue and the None in grey. Data labels indicate the absolute amount of state reaching the attractors. None of initializations reached an alternative attractor, even for higher amount of random initializations The number of initialization has no significant impact on the basal canalization and 10.000 initializations were considered sufficient to screen the state space in the current study.

Additional file 7: Table S3. Validation of transcriptional data-inferred regulatory interactions. Transcriptional interactions were inferred from the merged OA dataset. Inference was run with three algorithms and only interactions that were present in the results of the three algorithms (GENIE3, ARACNE, TIGRESS) were considered as additions to the model. An interaction was considered present for one algorithm if it scored higher than a threshold defined as the difference between the mean and standard deviation of all scores (Additional file 6, Data S2 ). Inferred interactions were validated whenever possible by looking for binding sites of the source transcription factor in the enhancer region of the target gene with GeneCards. The complete list of GeneHancer IDs provided by GeneCards are reported in the last column. ‘NA’ indicate that no binding site was found. The exact gene names and gene IDs that were queried in GeneCards are reported in the first columns.