dataesr/scientific-paragraphs-categorization huggingface dataset:
Specifically, Brodtkorb and Saetra use multiple GPUs to solve the Euler equations, whereas Nordhagen et al. use CPU clusters for simulating interacting particles by variational Monte Carlo methods. Software heterogeneity accompanies hardware heterogeneity, because new hardware architectures typically require completely new programming strategies. For GPU computing, although CUDA is a mature programming model, the low-level programming that is required is still too demanding for many domain scientists. Lifting the level of GPU programming is thus desirable. Brodtkorb and Saetra present a high-level programming framework for implementing a Euler equation solver in the Python programming language. Except for the core numerical computation that is implemented in CUDA C++ (and compiled just-in-time), the entire code is written compactly in Python.
Specifically, new Python classes are developed for managing the subgrids and individual MPI processes (each controlling a GPU). Moreover, the existing Python modules PyCUDA and mpi4py further alleviate the programming challenges associated with parallelization. In the same spirit, Markidis uses high-level programming interfaces, in particular, TensorFlow Quantum  [8]  and Strawberry Fields  [9] , to abstract the specifics of quantum hardware, thus avoiding explicit management of data movement or offloading of the quantum circuit execution. At the other end of the programming spectrum, Burchard et al. are forced to use the low-level Poplar C++ library that is specific for Graphcore IPUs. In general, new research is needed to overcome the programming hurdles that are often associated with heterogeneous computing.
Cellular senescence is a cell fate that prominently impacts physiological and pathophysiological processes. Diverse cellular stresses induce it, and dramatic gene expression changes accompany it. However, determining the interactions comprising the gene regulatory network (GRN) governing senescence remains a challenge. Recent advances in signal processing techniques provide opportunities to reconstruct GRNs. Here, we describe a GRN controlling senescence integrating time-series transcriptome and transcription factor depletion datasets. We infer a set of differential equations modeling the CS transcriptome using the "Sparse Identification of Nonlinear Dynamics" (SINDy) algorithm, discriminate genes with potential hidden regulators, validating the inferred GRN for time points not included in the training data.
Further to that, SINDy can be applied to distinct data subsets to identify possible HSs. However, SINDy and HS analysis have not been systematically used for gene expression modeling despite their demonstrated potential for modeling. In this work, we integrate time-series experiments of cells undergoing oncogenic RAS-induced senescence (OIS)  (Martínez-Zamudio et al. , 2020)  to build a system of differential equations that describes the GRN governing OIS using SINDy and highperformance computing. We validate the GRN generalization by comparing simulation results with experimental data acquired after the inhibition of two regulators of the CS phenotype, AP1-cJUN, and RELA, showing that model performance is dependent on the functions performed by each gene. Finally, we highlight the GRN genes whose dynamics can be inferred by hidden source analysis.
To empower our modeling approach further, we included transcriptomic data sets from RAS-OIS cells knocked down (KD) for the expression of cardinal transcription factors (TF) AP1-cJUN and RELA at 72h and 144h and ETS1 at 144h after RAS-induction by siRNAs collecting samples 48h after KD as previously described (Supplemental Figure  S1A )  (Martıńez-Zamudio et al. , 2020). Specifically, we used 72h AP1-cJUN and RELA KD and 144h ETS1 KD data sets for model training and the remaining two 144h KD datasets for AP1-cJUN and RELA for validation. The inclusion of TF KD datasets provides two benefits: first, the "transcriptome landscape" covered by the model training data is expanded, thus increasing the range of predictions it can perform  (Casadiego et al. , 2017). Second, it increases model resolution.
From a data-driven perspective, genes with highly correlated expression profiles in one condition provide the same information to the model and cannot be individually distinguished. Therefore, by integrating disparate datasets, it is likely that the expression of highly correlated genes develops distinct responses, allowing the model to identify covariates that regulate each gene separately and refining the inferred gene regulatory network (GRN). Together, we used these data sets for model training to compute model coefficients and model validation, as outlined below. Integration of disparate datasets naturally necessitates normalization to ensure that the effects being accounted for by the model are faithfully reflected in the biological process and not by technical variations inherent to high-throughput technologies.
Therefore, as a first step, we normalized all datasets for technical variation using ComBat  (Johnson et al. , 2007). Next, we independently clustered each time course using WGCNA 5 This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn. com/abstract=4199883 P r e p r i n t n o t p e e r r e v i e w e d  (Langfelder and Horvath, 2008) , identifying eleven to seventeen gene modules with highly correlated expression trajectories (Figure  1B  and Supplemental Figures  S1B-D ). 5124 genes were differentially expressed (> 30 % in expression levels) compared to un-induced, proliferating cells at 0h across the four experimental data sets. We aggregated gene sets belonging to the same module for all training datasets, i. e. , genes with expression profiles highly correlated in all experimental conditions.
We tabulated the WGCNA module for each gene in each experimental condition and identified seven hundred and ninety-nine distinct transcription module combinations (Figure  1C ). Each combination represents a unique temporal expression profile. We assigned each transcription module combination to a variable from which dynamics was inferred mathematically in this work. Three hundred seventy-five combinations are associated with one gene (Supplemental Figure  S1E ). Since AP1-cJUN, ETS1, and RELA are associated with variables corresponding to more than one gene, we assigned these TFs to an extra variable, resulting in eight hundred and two variables for our equation system. Altogether, we identified distinct gene expression profiles and associated variables using time-resolved transcriptomes from RAS-OIS cells depleted for TFs AP1-cJUN, RELA, and ETS1 that form the foundational basis for the inference of a dynamic mathematical system for RAS-OIS transcriptome modeling.
The latter is inferred independently for each variable by performing a sparse regression in the context of the CoS paradigm. We implemented the SINDy algorithm using the R statistical language (R Core Team, 2017) integrated with the Tensorflow library  (Abadi et al. , 2015)  and run in a High-Performance Computing facility containing Graphics Processing Unit (GPU) clusters. The produced scripts are available as R Markdown files and can be executed by a Singularity container, including all the required libraries available at ZENODO (see section Data and code availability). After identifying the coefficients for each of the eight hundred and two variables in our mathematical system (see Figure  1C ), we determined model performance by simulating the obtained equations for the RAS-OIS time-course data set, depicting in a heatmap observed and simulated time-resolved expression profiles (Figure  2B ) and in a corresponding histogram the correlation coefficients between expression profiles for each line  (Figures 2C).
This analysis demonstrated a high correlation between the experimental and simulated data with 429/802 (53 %) variables resulting in a correlation coefficient greater than 0. 9 (Figure  2C ). Next, we extended model performance to the AP1-cJUN (72h), RELA (72h), and ETS1 (144h) KD training datasets. In line with the RAS-OIS results, for the AP1-cJUN KD training set (Supplemental Figures  S2B, C ), we observed a significant congruency between experimental and simulation results, with 662/802 profiles (83 %) having a positive correlation coefficient greater than or equal to 0. 9. Similarly, for RELA KD (Supplemental Figures  S2D, E ), we obtained 583/802 profiles (73%) having a positive correlation coefficient greater than or equal to 0. 9. Finally, for ETS1 KD (Supplemental Figures  S2F, G ), we obtained 737/802 (92%) profiles having a positive correlation coefficient greater than or equal to 0.
9. To validate the inferred equations, we simulated the 144h AP1-cJUN KD (Figures  2D-E ) and RELA KD transcriptome datasets (Supplemental Figures  S2H-I ). Figure  2D  shows a heatmap for the experimental and simulated results of AP1-cJUN KD, their correlation for each variable, and the corresponding histogram with the correlation coefficient distribution (Figure  2E ). Most profiles (472/802=59%) have a positive correlation coefficient greater than or equal to 0. 9, while the minority (99/802=12%) have a negative correlation coefficient smaller than or equal to -0. 9. For RELA KD, 298/802 profiles (37%) exhibit a positive correlation coefficient greater than 0. 9 and 193/802 (24%) a negative correlation coefficient smaller than -0. 9 (S2H-I), suggesting RELA transcriptomic response might require more experimental data to be accurately predicted.
Hidden source analysis highlights genes with inferable dynamics Apart from allowing the reconstruction with a lower minimum on the required number of experimental time points, the CoS paradigm is also helpful in identifying nodes with hidden sources (HS), i. e. , regulators whose activity is not included in the input data  (Su et al. , 2014). We used this strategy to estimate variables with retrievable dynamics given our available datasets in this work. The hidden source (HS) analysis consists of using SINDy for distinct subsets of the time-course data and assessing the variance of each coefficient for data subsets at different time points  (Wang et al. , 2016; Shen et al. , 2014; Su et al. , 2014). Coefficients showing increased variance for a high proportion of the time-course suggest that HSs regulate the variables they describe and that the available data lack necessary information for their modeling.
To identify the biological functions of genes corresponding to variables associated with or without HS, we performed a pathway enrichment analysis using the Molecular Signatures Database (Figure  3B ). We observed that many pathways enriched in genes associated with variables without HS are also enriched in genes potentially associated with HS, including E2F target genes and genes downregulated in ultra-violet response. These pathways are also enriched when considering all differentially expressed genes in cells undergoing RAS-OIS with or without TF perturbation. We also find that genes corresponding to variables without HS contain a higher ratio of MYC targets when compared to all other genes. MYC is a master TF controlling a broad set of biological processes, including protein synthesis, metabolic homeostasis, transcription, cell proliferation, and tumorigenesis  (Chen et al.
, 2018; Dang, 2012). By contrast, genes related to variables with hidden regulators are enriched for pathways, including TNFα signaling via NFκB, epithelial-to-mesenchymal transition, mitotic spindle, and genes upregulated in KRAS signaling. To determine the accuracy of the inferred equations relative to genes without HS, i. e. , variables displaying decreasing σavg, we simulated the model for the AP1-cJUN and RELA KD transcriptome validation datasets at 144h. At each iteration of the equation solver, the next transcriptional state of genes without HS was computed by the equations inferred from the SINDy algorithm. The gene expression values of variables connected to HS were retrieved from the interpolated experimental datasets. The simulation results for initial conditions corresponding to AP1-cJUN inhibition 144h after RAS-OIS induction are shown in Figures  3C  and 3D.
Remarkably, the obtained equations predicted the effect of AP1-cJUN silencing at later time points in CS that were not included in the training dataset. At the same time, this was only partially the case for RELA silencing. As previously published (Martinez-Zamudio et al. , 2020), AP1-cJUN is a master regulator of the senescent phenotype, whereas RELA is a downstream effector TF of limited transcriptional impact. Because the simulation had a lower performance for RELA implies that RELA-mediated transcript output depends on hidden factors not only represented in the transcriptomic layer. The latter also underscores the need to accurately leverage diverse experimental data sets to describe TF function and transcriptional impact accurately. Given that TF knock-down data came from only twotime points, most profiles are monotonic, depending on the direction of the interpolation curve on the day of treatment.
The latter explains why the correlations for the TF depletion experiments had high amplitudes and were either positive or negative. To evaluate the influence of transcription regulatory events not present in our datasets, we applied a hidden source (HS) analysis  (Wang et al. , 2016; Shen et al. , 2014; Su et al. , 2014). This procedure measures the variance of coefficients when portions of data with distinct sizes are used to infer them. Thus, variables displaying an increasing variance with more data are subject to regulators not represented in the dataset. In the scope of this study, those hidden regulators might correspond to transcriptional regulators including post-translational TF modifications, histones modifications, or even co-factors of epigenetic modifiers  (Mitsis et al. , 2020).
Limitations of the study To compute the expression rate at each time-point accurately, we interpolated experimental data. This could be a potential source for artifacts for genes with rapidly changing time profiles. To address this issue, we will infer expression rates from single-cell RNA sequencing data (La  Manno et al. , 2018)  in the next iteration of the model. A sufficiently comprehensive dataset will allow us to define one specific time profile for each gene, eliminating the need to map multiple genes to a single model variable and increasing model resolution even further than our current solution. Another limitation of our model is that it scores only differentially expressed genes in at least one sample in the training dataset. By combining distinct experimental conditions, including assays where potential genes of interest present at least one activated and one inactivated observation, our approach will predict the impact of their perturbation.
STAR Methods Microarray transcriptome data download, preprocessing, statistical analysis, and annotation. The raw Affymetrix HTA 2. 0 data was downloaded from the Gene Expression Omnibus database (BioProject PRJNA439263, accession numbers GSE112084 and GSE143248;  Martínez-Zamudio et al. , 2020)  and pre-processed using Bioconductor R packages. All samples were normalized using the robust multichip average (RMA) tool implemented by the oligo R package  (Carvalho and Irizarry, 2010) , and batch effects were removed using the Combat tool  (Johnson et al. , 2007). Affymetrix probes were annotated using the hta20sttranscriptcluster. db R package (MacDonald, 2017) and internal control probes were removed. Genes with differential expression lower than 30 % compared to uninduced cells were filtered.
Hierarchical clustering and identification of unique expression time profiles The training datasets, consisting of the Ras-induced senescence time course (RAS OIS), the inhibition of AP1-cJUN and RELA at 72h, and the inhibition of ETS1 at 144h after RAS activation, were aggregated as shown in Figure  1A. The genes in each dataset were clustered independently with the WGCNA tool  (Langfelder and Horvath, 2008) , where the median of its replicates represented each sample. The parameters minimum cluster size, deepSplit, and threshold for merging clusters were set to 100, 3, and 0. 85. The soft threshold parameter was determined for each dataset separately, with the choice of the lowest value leading to a high Scale-free topology fit by applying the elbow method. Genes sharing the same modules in all time courses were aggregated, constituting one variable in the model.
Model inference 17 This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn. com/abstract=4199883 P r e p r i n t n o t p e e r r e v i e w e d To satisfy the RIP, the time course samples were scaled to zero mean and unit variance  (Wang et al. , 2016; Candes et al. , 2006). To compute the derivative of each variable time profile, each time course was smoothed by interpolating it by a factor of 24, and the respective derivatives were numerically computed using the secant method. The time profiles for each dataset were concatenated into a single matrix X, with rows representing the normalized expression levels for a given time-point and columns defining one variable in the model (Figure  S2A ). The time derivatives for each dataset were also concatenated in a single array ẋ i per variable i, such that each element matches the corresponding row in X.
Model simulation All model simulations were run with the R package deSolve  (Soetaert et al. , 2010). At each solver iteration, the model's current state was concatenated with another array consisting of its own pairwise multiplied elements, analogous to the computation of the Θ(X) matrix described in the previous section. This expanded array was multiplied by the inferred Ξ matrix, yielding the variation in gene expression given a specific state.
Hidden sources analysis To assess how many and which genes present a dynamics that can be reliably inferred from the available data, we performed a hidden source (HS) analysis as described by  Wang et al. (2016) ,  Shen et al. (2014), and Su et al. (2014). The matrix Θ(X), containing 284 rows, was randomly sampled five times for four distinct data points ratio Rm (20 %, 40 %, 60 %, and 80 %), and the inference procedure described in section "Model inference" was run independently for each sub-sampled matrix. For each value of Rm, the coefficient variance was computed, and the average variance (σavg) for all coefficients describing the dynamics of a single variable was assessed. Variables with decreasing σavg with respect to Rm were considered to be regulated by genes with expression levels accessible in our datasets and therefore, be described by deducible equations  (Wang et al.
, 2016; Shen et al. , 2014; Su et al. , 2014). To ensure a strictly decreasing σavg, we only considered variables where σavg for a specific Rm value is 80 % lower than σavg for the immediate lower Rm value. To validate the performance of the predicted model when considering only the variables without HS, we simulated the model with initial conditions corresponding to the validation datasets, i. e. , the KD of AP1-cJUN and RELA in senescent cells at 144h after RAS induction. The simulations were performed as described in section "Model simulation," where the expression levels corresponding to variables with non-decreasing σavg were retrieved from the experimental time courses. The performance of each simulation was assessed as described in section "Model assessment. ".
Pathway enrichment analysis Aiming to identify signaling pathways enriched for genes associated with variables with and without HSs, we performed an over-representation analysis for both groups using the msigdb  (Dolgalev, 2020)  and clusterProfiler  (Wu et al. , 2021 ) R packages, combined with the Molecular Signature Database (MsigDB) hallmark gene sets  (Subramanian et al. , 2005; Liberzon et al. , 2015). This analysis was also repeated for the set of all genes with at least a 30 % change in expression during all time courses.
Inference tools benchmark To evaluate the performance of distinct optimization methods on different platforms, we used the example provided by the R1magic R package  (Suzen, 2015). An array y is computed from the product of a Gaussian Matrix Φ and a sparse array x0. The execution time for each tool to recover x0 was registered for five distinct combinations of Φ and y. The evaluated optimization implementations were (a) the L1 regularization implemented by the R1magic R package  (Suzen, 2015) , (b) a BFGS implementation for 20 This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn. com/abstract=4199883 P r e p r i n t n o t p e e r r e v i e w e d CPUs (R Core Team, 2017), and GPU-compatible implementations of the (c) FTRL  (Abadi et al. , 2015)  and (d) BFGS algorithms  (Virtanen et al.
Data and code availability All data produced and software implemented are available on Zenodo (https://doi. org/10. 5281/zenodo. 6521071). Available data includes every file required to reproduce the results in the manuscript, including the generated time-course matrices; their nonlinear expansions; model coefficients considering the whole dataset and subsets used for HS analysis; and simulation results. The software implemented is available as both Rmarkdown files and html reports. We also share a Singularity container that runs model inference in a High-Performance Computing facility. Previously published transcriptome data are hosted in the GEO website under accession codes GSE112084 and GSE143248. 22.
Supplemental information titles and legends This preprint research paper has not been peer reviewed. Electronic copy available at: https://ssrn. com/abstract=4199883 P r e p r i n t n o t p e e r r e v i e w e d Table  S1 : Model variables metadata. For each DE gene identified, we provide metadata including microarray probe ID, Entrez ID, and gene symbol, in addition to the index of the variable in the model, the transcriptional module for each training dataset, the association of the variable to HS and the transcriptional module the same probe was associated to in our previous study  (Martinez-Zamudio et al. , 2020).
Model simulation with the training dataset. Heatmaps depict model-observed (left) and -simulated (right) time-resolved RAS-OIS gene expression profiles at indicated time points for each model variable. Correlation coefficients between observed and simulated gene expression profiles are shown as a bar-plot for each line in the graph at the far right of heatmaps. Data are expressed as row Z-scores. [C] Histogram depicting Pearson coefficient correlation distribution between model-observed and -simulated gene expression profiles of RAS-OIS transcriptome time-course for each model variable. [D] Model simulation with validation dataset. Comparison of observed (left) and -predicted (right) RAS-OIS gene expression profiles at indicated time-point for control and AP1-cJUN knock-down. Correlation coefficients between observed and simulated gene expression profiles are shown as a bar-plot for each line in the graph at the far right of heatmaps.
Acknowledgments J. A. N. L. F. d. F. was supported by  La Ligue Nationale Contre le Cancer ,  Association Nationale de la Recherche (ANR)  and  INSERM-AGEMED. O. B. was supported by  CNRS  and  INSERM-AGEMED. We thank past and present lab members, especially  Bertrand Friguet , for support and valuable discussions. The authors acknowledge support by the  High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen , the state of  Baden-Württemberg  through bwHPC and the  German Research Foundation (DFG)  through grant no  INST   37/935-1 FUGG.