dataesr/scientific-paragraphs-categorization huggingface dataset:
Implementation Designed as an R Bioconductor package, GWENA is a modular pipeline intended to ease the construction, interpretation and comparison of GCN. It reproduces a classical GCN analysis reinforced by complementary tools (Figure  1 ). ○ Co-expression network construction : computation through modified WGCNA function of a correlation matrix on the gene expression matrix, then transformation into an adjacency matrix, and finally into a topological overlap matrix (TOM). 4 ○ Modules detection : genes clusterization over the TOM with another modified WGCNA function. 5 ○ Biological integration : gene set enrichment of each module using g:Profiler services, and phenotypic association if a phenotype matrix is provided to describe the samples. 6 ○ Graph analysis : transformation of the TOM in a graph to compute different topological metrics and detect the hub genes.
-Filtering Genes are not always informative for modules detection as genes not always vary and their expression can be linked to technical biases. An additional filtering step can thus be applied to avoid noise and speed up the pipeline analysis. This operation must be carried out with caution as it may impact the network construction. Over-filtering may result in loss of informative signal and changing the data distribution could break the scale-free topology  [25, 23]. In addition, co-expression network analysis is a method designed to handle larger amount of data than differential expression analyses and capture the slightest significant genes expression variations. Two filters meeting these criteria are available in GWENA: 1 Low count filter : removes genes having a lower count than a pre-defined threshold (default is 5).
-Co-expression network construction The well-known R package WGCNA  [3]  was modified in order to be integrated it in our modular pipeline : the co-expression network construction step which computes the genes pairwise co-expression score has been isolated in its own R function. The first step of the co-expression score computation is the calculation of a correlation matrix based on the gene expression matrix. The Spearman correlation was added to the automated version of network construction in WGCNA as it ensures a better representation of genes monotonic relationships  [26]. A power law distribution is then fitted on the correlation matrix and the estimated power is used to elevate the correlation matrix, resulting in an adjacency matrix  [27]. According to the hierarchical organisation of gene co-expression networks  [28] , a topological overlap matrix (TOM)  [27]  is then computed using the adjacency matrix which represents the final gene co-expression score matrix.
-Biological integration Biological integration consists of two different analyses, namely gene set enrichment and phenotypic association. The gene set enrichment (or functional enrichment) analysis is performed using g:Profiler  [29]  through their gprofiler2 R package. Their enrichment function covers 9 biological functional databases: Gene Ontology (GO)  [9] , Kyoto Encyclopedia of Genes and Genomes (KEGG)  [30] , Reactome  [10] , Transfac  [31] , miRTar-Base  [32] , Human Protein Atlas (HPA)  [33] , CORUM  [34] , Human Phenotype ontology (HP)  [35] , WikiPathways  [36]. Realizing a custom enrichment file through a Gene Matrix Transposed (GMT) format in gprofiler2 requires the use of additional functions. Also, gpro-filer2 does not provide a merging function between the output of classical and custom enrichment to return all the enrichments in a single output.
-Modules differential co-expression Analysis of module preservation or unpreservation can be performed between different conditions such as treatments or phenotype. To isolate modules whose topology changes between conditions, GWENA first performs a permutation test using the NetRep R package  [38]. Seven topological metrics are computed on each module in each condition. A permutation is then applied to the selected control condition where each node label of the modules is randomly reassigned without replacement to another and the seven metrics are then recomputed on it. Using these permutations as a null distribution  [39] , modules are considered preserved if all seven topological metrics are significant for the alternative hypothesis (one or two sided) with the chosen alpha error. As the unpreservation of a module cannot be assumed from the non significant modules, a second step of preservation evaluation is carried out using a Z summary score  [40, 41].
Single condition modules analysis To illustrate the process of analyzing a single condition with GWENA, we first focused on studying the muscle gene co-expression computed in the young sub-population. The 95 modules detected on the coexpression score matrix with GWENA were merged according to their similarity indices, which resulted in a total of 35 modules (Figure  2. A. ). Each module was then tested for its association with a selected set of phenotypes related to muscle aging (i. e. age, sex, ethnicity, body weight and BMI) to isolate modules of interest. As shown in Figure  2. B, 15 of these modules were significantly associated with at least one of the phenotypes. These modules were provided to GWENA enrichment analysis (p value <0. 05 with g:SCS multiple testing correction) to identify their biological functions and assess their potential involvement in muscle function (Table  1 ).
All modules were at least enriched in one term and 8 obtained enrichment terms related to muscle activity or metabolism (Table  1 ). Modules 19, 21 and 25 were the top 3 enriched for terms related to muscle function. However, modules 21 and 25 terms were mostly coming from Human Protein Atlas and were also related to a wide range of additional tissues such as the pancreas, the cervix, the bladder, the stomach, or the skin and were thus deemed less specific for muscle aging than module 19. The number of enrichment terms is the sum of the terms across all sources of biological databases. The ratio of terms linked to muscle is the number of terms containing a term in a defined corpus of regex over the total number of enrichment found. The corpus is composed of the following regex: "muscle", "sarco\\w+", "myo\\w+", "muscul\\w+", "actin\\w+", "myosin\\w+".
Briefly, the remaining module 19 presented 77% of genes positively correlated to its eigengene (therefore 23 negatively, Figure  2. D. ), and the muscle enriched terms involved muscle adaptation and negative regulation of hypertrophy (Table  2 , Figure  2. C. ). The detection of hub genes by GWENA returned 12 hub genes, some of which are known as transcription factors. Among them, ARID5B (ENSG00000150347) is a transcription factor strongly co-expressed with KLF15 (ENSG00000163884) and TRIM63 (ENSG00000158022) (Figure  2. E). These two genes are present in the GO term GO:0014888 (striated muscle adaptation) to which ARID5B is not associated. The function of ARID58 is well known in adipocytes and hepatocytes but is still rarely studied in skeletal muscle metabolism. However, the knockout of this gene in mice has shown structural defects in the sarcomere structure  [45].
Multiple conditions modules comparison and analysis Differential expression analysis allowed to detect genes involved in aging in the last years (GenAge  [49] , Digital Ageing Atlas  [50] ). Such discriminant analysis is limited in helping to understand aging as this phenomenon is composed of concomitant mechanisms  [51]. Understanding the interaction between the genes is therefore crucial to determine the altered functions and the changes involved. Differential GCN between conditions overcomes this problem by detecting the subtle pattern modifications. Using our previously defined young (20 to 30 years old) and old (50 to 60 years old) skeletal muscle modules, we ran GWENA's GCN differential co-expression functionality to compare the modules between these age ranges. The GCN of each module detected in the young sub-population were taken as a reference and tested against the ones detected in the old sub-population.
Genes in sub areas of the module's network in the younger age range would then disconnect from some genes to reconnect to other genes in the old age range, forming different sub areas. Going beyond GWENA's functionalities, we investigated this reorganization of the gene connections occurring between the young and old network in the module 7. The reorganizations in sub areas were investigated through a partitioning around medoids (PAM) clustering method  [59, 60]  on the filtered co-expression matrix (0. 95 score filter). The clustering applied on both conditions returned 5 optimal sub-modules for the young, and 6 optimal sub-modules for the old. The gene composition of the sub-modules was highly similar between the condition (at least 81% common genes). Most of the differences in the gene composition are due to the disconnection of peripheral genes as previously spotted, and a small portion of the differences are due to the reconnection of genes or their attribution to another sub-module (Figure  4.
A and B). The new sub-module (sub-module 6, 4. B) appearing in the old age range is especially interesting as its creation is at the expense of the sub-module 1 of the young age range. Of the 13 genes composing it, 8 of the genes are from the sub-module 1, 3 are reconnecting genes, and 2 are from sub-modules 2 and 3. A gene set enrichment analysis with GWENA of this submodule 6 revealed significant enrichment in functions related to wound healing, coagulation, vessel diameter, platelet degranulation, and plasminogen activation (Additional file 1. 5). This is coherent with known morphological alterations of the vascular system in the aged skeletal muscle  [61, 62] , and the global immune/inflammatory response increased in aging  [54]. Also, 51 enrichments from young sub-module 1 were not found significant in any of the sub-modules in the old age range.
These enrichments involve antibacterial humoral response and negative regulation of endopeptidase activity. These terms are known to be associated with satellite cells (muscle stem cells responsible for muscle regeneration) regulators released by the vasculature in higher quantity in young skeletal muscle  [62]. To complement these analyses, we investigated the variations of co-expression scores leading to the appearance of sub-module 6 in the old age range. Using the network co-expression matrix δ returned by GWENA for each condition, a co-expression difference matrix (Figure  4. C) was computed such as δ oldδ young. In this matrix, gene pairs with a negative score indicates a decrease in the co-expression over aging while a positive score indicates an increase in the coexpression. Among the variations, 3 genes showed a significant increase in co-expression between them but also towards other genes of sub-module 6.
GWENA's contribution and comparison with existing tools Weighted GCN can be computed from existing tools such as WGCNA  [3] , wTO  [65] , CEMiTool  [66]. As both GWENA and CEMiTool use elements from WGCNA, they share notable functionalities. They use similar network construction and modules detection functions from WGCNA but offer their own filter on the datasets. On its side, wTO used a different version of a topological score to construct the network as they don't perform a power law conversion on the correlation matrix and don't use the same definition of topological transformation. Therefore, the main differences between WGCNA, wTO, CEMiTool and GWENA lie in the added functionalities for module analysis. Regarding biological integration, wTO provides neither phenotypic association nor gene set enrichment (Table  4 ).
The other three tools allow phenotypic association but differ on gene set enrichment analysis. While CEMiTool only allows enrichment on imported GMTs, WGCNA and GWENA allow enrichment on gene ontology. GWENA also allows enrichment on other databases of pathways, regulatory agents, and proteins (in addition to imported GMTs). Additional topological analysis functions are also available in several of these tools. The most common, hub gene detection, is present in WGCNA, CEMi-Tool, and GWENA in different forms. CEMitool and WGCNA offer respectively as hub gene the top 10 most connected genes and genes with a top kME score (membership module based on eigengene). However, methods based on a fixed number of hub genes tend to bias the information since the number of hub genes can vary according to the number of genes present in the module.
GWENA therefore proposes several methods (highest connectivity, superior degree, Kleinberg's score) based on a selection of genes with a hub score above a threshold. Finally, GWENA includes a differential co-expression analysis in the analysis pipeline as opposed to packages dedicated solely to it (DiffCoEx  [67] , CoDiNA  [68] , CoXpress  [69] ). The method in GWENA differs from the one present in WGCNA in that it includes a permutation test to prevent the problem of multi-testing. With the addition of the Z-summary score to detect unpreserved modules, GWENA is therefore the only pipeline including a differential co-expression analysis with high confidence in modules found unpreserved. GWENA, as other GCN analysis tool has limitations. A first one common to all GCN construction method is that the quality of input data (e.
g. filtration and/or proper normalization) will inevitably bias the results, especially if it breaks the scale-free property. A second limitation is the design of the permutation test that prevents reporting a significant unpreservation. The non-rejection of the null hypothesis of unpreservation can only state a lack of evidence of preservation. Therefore, unpreserved modules are determined among these modules lacking evidence (the non-significant modules) by the calculation of Z summary which only provide a tendency in the unpreservation  [40]. The present application of GWENA to skeletal muscle aging also presents its own limitation. All analyses were performed on skeletal muscle sample and results were commented regarding this context. However, to be sure of the specificity of the findings, an additional differential co-expression of the modules should be performed on samples from other tissues from subjects with similar age range.
Conclusion In this paper, we introduced GWENA, an R package on Bioconductor to construct and analyze GCN in a single pipeline through a whole range of tools from biological integration, topological analysis, and differential co-expression. The package eases the GCN analysis through simple input and output functions combined to a set of visualizations to explore the results. The separation of each step of the analysis in one function also allows quick and easy replacement if users wish to use another method for this block. GWENA demonstrated its performances on both single and multiple condition analysis through an exploration of variations of skeletal muscle function and processes in aging. The single condition analysis showed it is possible to find new genes potentially involved in an existing GO annotation using hub genes, network neighboring genes and enrichments.
Consent for publication Not applicable Figures Figure 1 Detailed steps of analysis performed in GWENA's pipeline, from expression data to characterization of the modules and comparison of conditions. 1 Input : expression matrix pre-normalized and aggregated to gene level if it is a transcript matrix. 2 Filtering : optional genes filtration according to transcriptomic input technology. 3 Co-expression network construction : computation through modified WGCNA function of a correlation matrix on the gene expression matrix, then transformation into an adjacency matrix, and finally into a topological overlap matrix (TOM). 4 Modules detection : genes clusterization over the TOM with another modified WGCNA function. 5 Biological integration : gene set enrichment of each module using g:Profiler services, and phenotypic association if a phenotype matrix is provided to describe the samples.
Availability of data and materials • Project name : Gene Whole co-Expression Network Analysis (GWENA) • Operating system(s): Platform independent • Programming language: R • Package on Bioconductor : https://bioconductor. org/packages/devel/bioc/html/GWENA. html. • Package GitHub repository (development, issues, and pull requests) : https://github. com/Kumquatum/GWENA. • Licence : GPL-3 • Source code for GTEx analysis : https://github. com/Kumquatum/GTEx v8 GWENA/. • GTEx v8 public data (RNA-seq and anonymized phenotype) : https://gtexportal. org/home/datasets. • GTEx detailed phenotype requires a dbGaP request.
In extreme cases, the FIM can even make some practically unidentifiable parameters appear as identifiable  [6]. As a consequence, the FIM is inadequate to study practical identifiability. Since our focus is on the practical identifiability of MEM, we will not use the FIM to assess the identifiability of our models. Other methods have been developed for models without mixed effects, such as the profile likelihood  [6] , but to our knowledge they have not yet been implemented in any of the existing software for MEM calibration. Given the widespread use of these software for the calibration of MEM in pharmacology and personalized medicine  [19] , it would be particularly interesting to be able to empirically study the practical identifiability of MEM, directly from the calibration software. Secondly, one might run the estimation algorithm several times, using a sample of initial guesses for the parameter values  [13].
Thus, the size of the sampled dataset (in terms of the number of individuals) is critical for the estimation of the population variances, since a sampling bias in a small dataset could mask the variance on some parameters. From the point of view of experimental design, the determination of the necessary sample size in order to guarantee a certain level of confidence on all parameters in MEM is a central question, which has already been covered to a certain extent  [22, 27] , using geometric features of the likelihood surface approximated from the FIM [28, section 10. 5. 3]. From the point of view of model reduction, it would be tempting to remove the unidentifiable variances by setting them to zero, potentially improving the estimation of the other parameters of the MEM without affecting the quality of the model fit to the data.
In some cases however, adding a random effect or a covariate to a MEM could improve parameter identifiability [26, Section 5. 1], as it might split out combinations of structurally unidentifiable parameters. In this paper, we adress these questions using a MEM of the in vitro erythropoeisis that we adapt from a previous model, proven to relevantly reproduce the dynamics of single replicates of the experiment  [29]. This MEM accounts for experimental variability by assigning different parameter values for proliferation and differentiation in each replicate of an identical experiment. We assess its identifiability using a multistart approach, based on extensive parameter estimations with the MEM calibration software Monolix  [30]. Then, we reduce the model in order to make it identifiable, using the correlations between the estimated parameter values.
The proportional error model introduces one additional parameter b i for each observable, resulting in three error parameters for our SCB model. Together with the dynamic model of System (1) and the parameter model of Eq. (  2 ), this error model defines our first version of a MEM for the in vitro erythropoiesis. Since all other MEM in this manuscript will have the same dynamic and error components, we will omit them from now on, and will define each MEM by its parameter equation only, such as Eq. (  2 ).
Indeed, our dataset only comprises 7 individuals. In order to obtain an identifiable model based on Model  (7) , one might remove the random effect on one or several individual parameters. Fixing their values across the population might allow for a more precise estimation of the other, still variable, individual parameters, while keeping the same quantitative fit as with Models (  2 ) and  (7). However, all parameters are not necessarily equivalent in this regard, since different parameter sensitivities would make the model output more flexible under some combinations of fixed parameters. This would allow for these combinations to better fit the data, depending on the sensitivies of the model output to the parameter values. This sensitivity is imposed by the analytical solution to the structural equations of the model that we defined in System (1), but for most models there is no closed-form expression for the parameter sensitivities.
However, the number of convergent SAEM runs is critical to the assessment of population parameters identifiability. Depending on the complexity of the model and the dataset, 50 SAEM runs might not be sufficient to assess parameter identifiability and allow for model reduction, and the number of runs to be performed should thus be finely tuned in order to avoid issues with computational time. Morevover, we used the correlations between population parameters to define Model  (7) , with constraints on the fixed effects. The exact shape of the likelihood landscape and the resulting unidentifiability is related to the structure of the model, and the quality of the data. This means that we were able to explore the parameter space near the likelihood optimum using pairwise correlations (Fig. 4 ).
Acknowledgements We thank the members of the  Systems Biology of Decision Making team at the Laboratory of Biology and Modelling of the Cell , as well as the  Dracula team of Inria , for insightful discussion all along the course of this project. We thank the  BioSyL Federation  and the  LabEx Ecofect  ( ANR-11-LABX-0048 ) of the  University of Lyon  for inspiring scientific events. We also would like to thank members of the Population Approach Group Europe that were present at the PAGE meeting 2019 for enlightening feedback. Finally, we want to thank  Pr. Saccomani  ( University of Padova, Italy ) for providing the manuscript for reference [9].
Geographical structure We use the freely available OSM database to extract all the points of interest for this study OSM. The accuracy of this database is partially high in France because the official national land registry (the 'Cadastre') database is merged into the OSM database. These databases are formatted in ASCII XML, with a size of approximately 80 GB for France. OSM labels residential buildings specifically. We use the surface of these buildings and their geographical position to allocate each agent to a 'home' (residency) building. For the mobility model, we also include all the other types of buildings (e. g. hospitals, schools, airports, commercial centres, etc. ), where agents can meet. In our current simulations, the initial database size contains 4. 1 10 8 nodes and 4. 8 10 7 buildings that need to be pre-processed to compute additional characteristics such as building surface, usage, or geographical location, which can be reused in different simulations.
Demographic model The Institut National de la Statistique et des Etudes Economiques (INSEE) publishes data with a high-level resolution of the distribution of the French population  (INSEE, 2019). Here, we use the information about the full population in each city in 2016 to allocate agents to different locations. By combining this database and OSM's, we can compute the number of residents in each building of each city. More precisely, agents are allocated to buildings proportionally to their floor-surface projection. The number of agents in each building is given by the equation N k ¼ N city S k P i S i þ a : (1) where N k is the number of agent in the building k, N city the number of residents in the city considered, S k the floor-surface projection of building k, P i S i the total surface of all residential buildings in the city, j Á j the entire part of a number, and a 2 ½0; 1 a scaling parameter such that P k N k ¼ N city.