dataesr/scientific-paragraphs-categorization huggingface dataset:
Data Availability All summary statistics are available at https://pdgenetics. org/resources. Four summary statistics tables are made available: (1) male PD GWAS (all data), (2) female PD GWAS (all data), (3) male PD GWAS (no UKB proxy data), and (4) female PD GWAS (no UKB proxy data). All code used has been made available on GitHub: https://github. com/neurogenetics/Autosomal-sexdifferences-PDv2.
Three different datasets were generated to feed our quantitative analysis pipeline (see Figure  1 ). The first one was used to estimate leaf initiation time for each leaf rank (Dataset 1). Plants were dissected and observed using a binocular microscope (Nikon SMZ645) to determine the number of leaves, including the smallest visible primordium. This ensured that leaves were counted as soon as they appeared as primordia at the periphery of the shoot apical meristem. For each genotype, the output was a collection of pairs (a p ,n) i , with a p the plant age when dissected and n the number of observed leaves for any individual plant i. Leaf number was determined in 295, 33 and 38 plants for Col-0, clfsep and kluh, respectively. The second dataset was used to estimate the temporal dynamics of leaf size (Dataset 2).
The third dataset was used for quantitative image analysis and reconstruction of growth trajectory in MorphoLeaf software (Dataset 3). For this, wild-type data previously generated were reused and enriched following the same acquisition protocol as previously  (Biot et al. , 2016). Using MorphoLeaf, following automatic leaf contour segmentation, the proximal limit of the blade contours were manually delineated and morphological landmarks were automatically computed  (Biot et al. , 2016). Manual corrections were applied if necessary. This yielded five sets of 2D images for Col-0 L3, L5, L7,  L9 and L11, with 160, 196, 168, 162  and 312 images, respectively, with organ length ranging from several tens of microns up to few centimetres for mature leaves. Images of L11 were similarly generated for the mutants, providing 214 and 390 images for kluh and clfsep, respectively.
Growth curves A method was setup to estimate the function relating organ size and age. Let us consider a given leaf in the rosette (e. g. , L11 in Col-0), for which initiation time T was previously determined (see Section 3), and observed pairs (a p ,l) i associating plant age and blade lengths (Dataset 2). Initiation time T was used to translate these data in time, to obtain pairs (a,l) i , where a = a p -T is the age of the leaf since its appearance in plant i. For any organ considered in this study, length evolution with time exhibited a sigmoid shape. The Hill model was the most efficient in the family of sigmoid functions to fit size measures (Supplementary Figure  S5 ). The function L relating length l and age a of the leaf from initiation was thus given by L(a) = L ∞ a n t 50 n + a n , (1) with L ∞ the length of the mature leaf, t 50 the time at which the length is half of L ∞ , and n the Hill parameter.
Morphogenesis quantification and reconstruction MorphoLeaf software was used to quantify leaf blade shape evolution during growth, from collections of 2D organ images (Dataset 3). Using segmented leaf contours and morphological landmarks, it provides measures of both global shapes (e. g. , leaf blade area or elongation) and serrations at the margin (e. g. , tooth area or aspect-ratio). The application also generates a mean growth trajectory, that is, a sequence of organ contours computed by timeweighted averaging of individual contours  (Biot et al. , 2016). Blade length was used as a proxy of leaf age in the initial MorphoLeaf release  (Biot et al. , 2016). Here, we modified the application to integrate real organ age. Using the estimated parameters of the growth function (L ∞ , n and t 50 ), MorphoLeaf software computes the age of each individual organ from its length [equation (2)].
Numerical and statistical methods Systems of differential equations were simulated within the COPASI software  (Hoops et al. , 2006)  using the RADAU5 integration method with default initial step-size and tolerance parameters. Model parameters were estimated in COPASI using the Evolutionary Strategy with Stochastic Ranking method  (Runarsson and Yao, 2000)  with default parameters. The cost function minimised in the estimation procedure was the sum of squared differences between predicted (N 0 ,. ,N K ) and experimentally obtained (M 0 ,. ,M K ) numbers of individuals with 0,. ,K organs at any time in the population. The experimental dynamics of the numbers M 0 (t),. ,M K (t) of individuals having 0,. ,K organs were not directly observed. Indeed, our experimental data (e. g. , Dataset 1 for the number of leaves per plant) were sets of couples t i ,k i , where t i was the age of the ith individual and k i its number of organs.
Temporal dynamics of leaf initiation We used our mathematical model to estimate the initiation times of leaves of different ranks (organ apparition orders in the plant) in Col-0. Using experimental recordings of the number of leaves as a function of plant age (Figure  2a ; Dataset 1, see Section 2), we computed the experimental dynamics of the proportions of plants with varying numbers of leaves in a typical population (Figure  2b , dots). The same procedure was applied to mutant data (Supplementary Figure  S4 ). Leaves of odd ranks only were considered in the three genotypes. In the wild type, we modelled the dynamics up to L11. In the mutants, we considered the dynamics from L9 to L11 (kluh) or from L7 to L11 (clfsep) to comply with the experimental samplings (Supplementary Figure  S4 A, B ).
Temporal dynamics of global leaf growth To translate blade length measures (Dataset 2) from plant age to leaf age, we first applied a temporal shift using estimated organ initiation time (Figure  2d ). The functions relating leaf size and age were then evaluated for wild-type and mutant leaves (Figure  2d -f and Supplementary Figure  S6 ). The resulting estimated functions revealed a strongly patterned, regular arrangement of leaf growth kinetics according to organ rank in Col-0 (Figure  2e ). There was a switch during growth in the ordering of leaf sizes of different ranks, while the final blade length increased with leaf rank, the reversed ordering was observed in the initial phase of growth. In addition, the graphs suggested that both the maximal growth rate and the time where it was reached also increased with leaf rank.
Leaf morphogenesis trajectories We next examined the morphological changes that accompanied growth dynamics. Following the delineation of leaf contours and extraction of geometrical landmarks on leaves collected at different timepoints (Dataset 3, see Section 2), a continuous morphological trajectory was reconstructed for each leaf rank using a new version of the MorphoLeaf software (Figure  3a ). The age since initiation of all collected leaves was computed using automatically measured blade lengths and growth kinetics functions. In a previous work, we compared developmental trajectories between leaves of different ranks by analysing stages of the same length  (Biot et al. , 2016). This was informative but growth kinetics functions (Figure  2e ) showed this induces comparing leaves at different ages.
Dating teeth To finely dissect the changes occurring during leaf morphogenesis, we next focused on individual teeth. To determine teeth initiation times, we applied the same strategy as the one used above for leaves. The number of teeth per half-leaf as a function of leaf age (Figure  4a , and Supplementary Figures  S9  and S10 ) was automatically determined with MorphoLeaf from segmented leaf contours of Dataset 3. We then fitted our mathematical model of initiation dynamics to the empirically derived proportions of leaves with different numbers of teeth. The model fitted closely to estimated teeth proportions at all leaf ranks and all tooth ranks (Figure  4b , and Supplementary Figures  S9  and S10 ). In some cases, there was some discrepancy between model and data for later teeth, probably due to some teeth missed in larger leaves, as they tended to smooth out  (Biot et al.
To model the epidemic spread through a trade network we use a stochastic compartmental model that takes into account demographic dynamics and animal exchanges. The population structure of the model is calibrated using real data. Our decision model is inspired by previous studies  [13] [14] [15]  , in which the result of a decision regarding the voluntary adoption of a control measure for a human Epidemic model with demography in a metapopulation based on a trade network. For this work we place ourselves in the context of a hypothetical livestock infectious disease that is transmitted only through a contact network structure consisting in herds that exchange animals. This population structure is inspired by real data on animal movements, extracted from the French Cattle Identification Database (FCID).
Epidemic-decision related Mean cumulative intra-herd incidence rate by aggregated vaccination pattern Vector of the mean cumulative intra-herd incidence rate (see output p [0,T] Ianim ) of herds grouped by the aggregated vaccination pattern: herds that never vaccinate, herds that vaccinate at least once and at most half of the time, and herds that vaccinate more than half of the time but not always. Without taking the first decision into account design  34. For all the sensitivity analyses we use the R package multisensi  35. In the PCAs the means are centered and scaled, and the dimension is selected as the smallest value that keeps at least 95% of the total variability. Among the many experiments, we retain the results of the three following ones: (i) First experiment: all 13 inputs. The means and variances of all outputs: by group, and all outputs simultaneously.
Acknowledgements This work was carried out with the financial support of the  French Research Agency (ANR) , project  ANR-16-CE32-0007  ( CADENCE). We are grateful to the  INRAE MIGALE bioinformatics facility (MIGALE, INRAE, 2020. Migale bioinformatics Facility , doi:  10. 15454/1. 5572390655343293E12 )for providing computing resources. We thank the  French minister of agriculture  for providing the access to the FCID.
Datasets compilations Most of the data collected are from the French Research Institute for Marine Exploitation's (IFREMER) REPHY (Monitoring network for phytoplankton and phycotoxins), RHLN (Hydrologic Network for the Normandy Littoral) and SRN (Regional Nutrient Monitoring) programs. The flora dataset gathers information on phytoplankton abundance and diversity on the sub-surface level at a bi-monthly or monthly resolution for the time period of  2008 -2019 (Neaud-Masson, 2016). The abiotic dataset either is from field sampling, or has been extracted from local weather networks (Météo France, Hydrological bank) on a daily A. Lefran et al. resolution.
Phytoplankton datasets The flora dataset was initially composed of 131 taxa (60% diatoms, 29% dinoflagellates) most of which were identified to the species or genus level or regrouped as sets of species or genus (85% of the taxa). In order to ensure robust analyses, taxa that represented more than 99% of the total abundance over the studied period and areas were selected. Moreover, only the abundances above 100 cells. L -1 were considered since this concentration is the detection limit with microscopy observations  (Neaud-Masson, 2016). Special attention was put towards the naming convention so as to ensure temporal coherence. The final list for the 108 considered taxa is given in the supplementary data (Table  S3 ).
Numerical analyses Data analyses were managed through R, version 3. 5. 1 (R Core Team, 2018), with the use of "ggplot2" package  (Wickham, 2016)  for the majority of the graphical outputs. Datasets and scripts are available on a GitHub repository (https://github. com/AngieLef/PhytoComm_Ch annel21). All phytoplankton cell counts were transformed to a logarithm base-10, with two decimals. Seasons are defined as three successive months with winter being the combination of the 12th month of year Y and the 1st and 2nd months of year Y+1.
Seasonal temporal trajectories of phytoplankton Trajectory analysis is an ordination that aims at comparing the behaviour of each community (trajectories) according to two main criteria; the resistance and resilience of communities, and the convergence or divergence for each year. For this, a recent Community Trajectory Analysis (CTA) was applied using the dissimilarity matrix described above  (De Cáceres et al. , 2019). It display, through annual arrows, the temporal evolution of the community assemblage over the study period. The function is part of the "vegclust" package  (De Caceres et al. , 2010)  and the analysis was made independently for each season (ie: all years included).
Linking environmental and community datasets can draw out the structuration of each stations' community, depending on the main abiotic influences. For this analysis, both data tables needed to be scaled to the same temporal index, namely for the sites and years without any gaps. This is why only the seasonal scale was considered. However, STCA station was missing a value for the N/P variable during Winter of 2019; therefore, the mean value of all past winters' N/P ratio was used to complete the dataset.
Acknowledgement We would like to thank the current and former Ifremer teams that have worked on the acquisition of this data and still are. Especially, thank you  Sylvaine Françoise  (Ifremer -Port en Bessin),  Aurélie Legendre  ( Ifremer-Dinard ) and  Pascale Hebert  ( Ifremer -Boulogne sur Mer ) for your help checking the naming convention. Also, the Quadrige 2 , REPHY, RHLN and SRN teams and networks that make the data easily available through their platforms and protocols. In addition, special thanks to all public data collectors within Méteo-France that helped us extend our data comprehension with abiotic information. Thanks to  Thibault Cariou  for helping with the map conception and  Carine Sauger  for the redaction review.
The alveolar macrophage host model  (Bordbar et al, 2010)  used to illustrate the metabolic interactions between multiple organisms was constructed from the Homo sapiens Recon 1 model  (Duarte et al, 2007)  based on previously built algorithms. This was followed by a set of manual curation processes to construct the context-specific alveolar macrophage host model  (Bordbar et al, 2010). Similarly,  Wang et al (2012)  developed several draft context-specific models using the Human Recon 1 model. We have built on these approaches by developing an integrated epithelial cell/SARS-CoV-2 metabolic model and used a combination of structural and dynamical analyses to assess the model and make predictions. We have applied the recently released community standards to facilitate our development of a standardised model for the systems biology international community and used the MEMOTE quality control software to assess and compare our model with previously developed GEMS  (Carey et al, 2020 , Lieven et al, 2020).
Construction of the human airway epithelial cell and the SARS-CoV-2 virus We constructed an integrated genome-scale metabolic model (GEM) of the human airway epithelial cell with the SARS-CoV-2 virus using the methods described in  Aller et al (2018)  and  Renz et al (2020). We used a previous draft reconstruction of the epithelial cell  (Wang et al, 2012)  and improved on the context-specificity of the GEM with the addition and simulation of data from the gene expression omnibus on COBRApy  (Ebrahim et al, 2013). We performed further refinement from literature and by using the human metabolic networks in the HumanCyc database  (Trupp et al, 2010)  and Recon3D  (Brunk et al, 2018). Model checks were done with MEMOTE, and refinements were performed with MetaNetX 4. 2, BiGG, ChEBI, MetaCyc, and PubChem databases  (Caspi et al, 2014; Hastings et al, 2016; Norsigian et al, 2020; Kim et al, 2021; Moretti et al, 2021).
The new GEM contains 4,660 genes, 3,614 reactions, and 4,052 metabolites and conforms to the minimum standardised content for a newly published GEM based on recently published community standards  (Carey et al, 2020) ; 100% of the metabolites in (iHsaEC21) have a human-readable descriptive name, 100% have an inchi key, 100% of metabolite annotation conformity with the BiGG database and in MetaNetX, Kyoto Encyclopedia of Genes and Genomes (KEGG), ChEBI, ModelSEED, HMDb, or MetaCyc  (Caspi et al, 2014; Hastings et al, 2016; Wishart et al, 2018; Norsigian et al, 2020; Kanehisa et al, 2021; Kim et al, 2021; Moretti et al, 2021; Seaver et al, 2021) ; 100% of the metabolites have a charge and chemical formula with a charge balance of 75. 3% (Supplemental Data 1). In addition, 97% of the reactions have a human-readable descriptive name, 100% of reactions conform with the BiGG database and as well as in MetaNetX, KEGG, ChEBI, ModelSEED, HMDb, or MetaCyc.
The enzymes of the reactions at the nodes of intersection are (1) host dependency factors identified by the model as necessary for the survival of the virus and (2) proteins from experimental datasets with high-confidence virus-human protein-protein interaction data; we refer to these points of intersection as PPi-pathway intersection nodes. The PPi-pathway intersection (P-Pi) nodes identified are present in biosynthesis pathways such as the cysteine and S-adenosyl-Lmethionine biosynthetic pathways. In both pathways, the enzyme S-adenosylmethionine synthase (MAT2B) catalyses the phosphorylation reaction of methionine to S-adenosyl-L-methionine. During infection, the viral protein Nsp9 is seen to react with MAT2B (Gordon et al, 2020) (Fig 3A  -C ). Another viral protein, Nsp8, also interacts with the enzyme O-phosphoseryl-tRNA(Sec) selenium transferase (SEPSECS), which catalyses the last step of the L-selenocysteine biosynthesis pathway (Fig  4 ).
After further refinements using Recon3D, HumanCyc, MetaNetX 4. 2, BiGG, ChEBI, MetaCyc, and PubChem databases  (Caspi et al, 2014; Hastings et al, 2016; Norsigian et al, 2020; Moretti et al, 2021)  including manual curations with literature, iHsaEC21 now has a MEMOTE quality score of 51% (Supplemental Data 1). Although the iAB-AMØ-1410  (Bordbar et al, 2010)  was also developed from the human Recon1 model, there was further refinement with literature, which attributed to its quality score of 45% (Supplemental Data 2) compared with the draft epithelial cell model from  Wang et al (2012)  which had a low-quality score of 19% (Supplemental Data 3). Previous studies have also demonstrated the role of the ACE2 as the receptor for both the SARS-CoV and the SARS-CoV-2. The ACE2 cells are expressed in the human airway epithelial cells  (Hamming et al, 2004; Liu et al, 2020).
GEM reconstruction, curation, and simulations We obtained an automated draft reconstruction of the human airway epithelial cell  (Wang et al, 2012)  and evaluated the metabolic functions and reactions collected in the draft reconstruction against organism-specific literature. We used the gene expression datasets of the human airway epithelial cell  (Deprez et al, 2019 Preprint; Vieira Braga et al, 2019)  to curate, annotate and improve on the model  (Fig 8). We manipulated and simulated the model, including the addition of the gene expression data, with COBRApy  (Ebrahim et al, 2013)  and GLPK (https://www. gnu. org/software/glpk/). We obtained additional reactions, gene-to-reaction associations, and pathways (that were not in the automated model) from HumanCyc  (Trupp et al, 2010)  to develop a new genome-scale metabolic model (GEM) of the human airway epithelial cell.
We improved on the new reconstruction by mapping the genes and reactions of the GEM to Recon3D, a standard human metabolic model, MetaNetX 4. 2, BiGG, ChEBI, MetaCyc, and PubChem databases  (Caspi et al, 2014; Hastings et al, 2016; Norsigian et al, 2020; Moretti et al, 2021)  and produced a revised human airway epithelial cell iHsaEC21 (Fig  8 ). The revised model reconstruction, iHsaEC21, can be instantiated without error on the COBRA software (version 0. 16. 0) (Ebrahim et al, 2013).
Model optimisation To improve on the quality of the model, we removed the reactions with dead-end metabolites, which were previously found in the automated model. This was followed by iteratively assessing the model and accounting for the cell-specific metabolic and exchange reactions in the epithelial cell model. We compared and performed various iterations of the new GEM with Recon3D, MetaCyc, MetaNetX 4. 2, BiGG, ChEBI, and PubChem databases  (Caspi et al, 2014; Hastings et al, 2016; Norsigian et al, 2020; Kim et al, 2021; Moretti et al, 2021). As a result, we added charge and formulae to 1,473 compounds on the new GEM. During the iterations and optimisation of the model, we performed quality control checks on MEMOTE, a standardised genomescale metabolic model testing software  (Leiven et al, 2020).
Integration of the human airway epithelial cell model and the SARS-CoV-2 virus We integrated the viral biomass maintenance function, previously developed for the macrophage cell  (Renz et al. 2020 (Renz et al. , 2021) )  into the new GEM to produce an integrated model of iHsaEC21+SARS-CoV-2. We interrogated the new model to identify the host-dependency factors for the SARS-CoV-2 virus by using the novel software tool findCPcli. The new GEM is encoded in the Systems Biology Markup Language (SBML) (Keating et al, 2020) (iHsaEC21); i for in silico, Hsa for H. sapiens, and EC for airway epithelial cell published in 2021. iHsaEC21 consists of 3,752 reactions, 3,914 metabolites, 4,660 genes, and 48 metabolic pathways. To assess and predict the performance of the models, we made use of FBA and FVA  (Orth et al, 2010).
For instance, if the growth rate predicted by FBA is μ max , then the range of fluxes of a given reaction i that are compatible with such growth rate can be where v growth is the flux of the reaction associated with growth and v i is the flux of reaction i. FBA and FVA were computed on the metabolic network of the host, both with and without the reaction modelling the production of the virus, by using the Python toolbox COBRApy. Bottleneck reactions, like chokepoint reactions, are required for the reaction synthesis and the removal of these reactions will cause an accumulation or depletion of the metabolites; thus, they represent potential drug targets. The software tool findCPcli was developed to compute bottleneck reactions on genome-scale models by considering the structural and dynamic information of the models.
Acknowledgements The authors of this manuscript are supported by The  Wellcome Trust   107032AIA  ( RA Floto, BP Bannerman ) The  UK Cystic Fibrosis Trust  ( Innovation Hub grant   001 ) (RA Floto), and the  Spanish Ministry of Science, Innovation, and Universities  ( J J úlvez  and  A Oarga ). The authors would like to thank  Jo Chukualim  for the graphics and illustrations. The authors would also like to thank  Charlotte Passemar  for slides with preliminary data demonstrating the expression of ACE2 in epithelial cells.
Data Availability The model has been deposited as an SBML Level 3 Version 1  (Hucka et al, 2018)  file with FBC extension and the minimal information required in the annotation of models (MIRIAM)  (Juty et al, 2012)  in BioModels (https:// www. ebi. ac. uk/biomodels/MODEL2007210001) (Malik-Sheriff et al, 2020). The source code of findCPcli together with its documentation is available at https://github. com/findCP/findCPcli. The tool requires Python 3. 5 (or higher) and can be installed with pip, the standard package installer for Python. The pathway maps were created with the pathway collage software  (Paley et al, 2016)  and the essential reactions are available as a smart table on HumanCyc (https://biocyc. org/group?id=biocyc17-29351-3833343490).