dataesr/scientific-paragraphs-categorization huggingface dataset:
It is possible to run the model with a few user data while relying on the underlying pedoclimatic, crop and fertiliser database. It is also possible to be much more precise by feeding the model detailed user data. The model should be used with system-specific data, and not with aggregated averages or generic data. A list of minimum data requirements is proposed in Table  6. Some of the data are available from the model implementation's database or generated after intermediate calculations, but results would gain in robustness if those would be personalised (highly variable and high impact on final results). equally between the crop time period receiving the losses and all the previous crops time period until the one with a drainage period". For a typical temperate field crop, this method allowed to share the nitrate losses during a crop time period equally between the crop and the previous crop.
Acknowledgements The authors are grateful to several colleagues which provided valuable discussion elements:  Montserrate Núñez ,  Erica Montemayor  and  Edilene Pereira  (IRTA);  Thierry Morvan and Virginie Parnaudeau (INRAE) ,  Laurent Thuriès  and  Samuel Legros  (CIRAD);  Julie Auberger  ( INRAE ) and the rest of the team behind the MEANS platform (https://pfmeans. inra. fr/means),as well as  Antoine Versini  (CIRAD),  Johan Blockeel  (FIBL) and  Marie Benoit  ( UPMC ) for the valuable datasets for model validation. The authors are grateful of the  Agroscope team in charge of SALCA-NO3 , for providing an Excel implementation, insight and valuable discussion:  Robin Harder ,  Maria Bystricky  and  Thomas Nemecek.
Software and data availability An operational implementation (prototype) in Excel of the Indigo-N v3 model, as well as all detailed model equations, are available from the Supplementary Material. All associated data is described in Galland et al. (submitted)  and available as a downloadable dataset: Galland, Victor; Avadí, Angel; Bockstaller, Christian; Paillat, Jean-Marie, 2020, "Data to inform the modelling of direct nitrogen field emissions from global agriculture", https://doi. org/10. 15454/RM6TYY, Portail Data INRAE, V1.
C, N, O abundances The abundances of C, N, and O are derived from CN, OH, and CO molecular lines. They are interdependent due to the molecular dissociative equilibrium. Since the molecular lines are spread all o v er the spectra, these abundances are derived first, and they are reported in Table  5. Computing synthetic spectra employing the PFANT code described in  Barbuy et al. ( 2018b ) , we have derived C, N, O abundances in two ways: a) In the region 15144-16896 Å, first we derive the O abundances by analysing the molecular lines of OH. Some of the most prominent OH lines in this region are at:  15264. 60, 15266. 160, 15278. 516, 15281. 045, 15719. 687, 15893. 524, 16074. 151, 16662. 187, 16872. 265, 16895. 164 (Piskunov et al. 1995 ; Ryabchikova et al. 2015 ) ,  Kurucz ( 1993 )  and NIST  (Martin et al.
Silicon and calcium Si and Ca are mainly produced during the e xplosiv e nucleosynthesis of SN II events (WW95; McWilliam 2016 ) with smaller contributions from supernovae type Ia (SNIa). The α-elements Si and Ca are plotted in Fig. 10  for the 58 sample stars together with literature data from Garc ía-P érez et al. (  2013   Barbuy ( 2017 ). Chemodynamical evolution models from Fria ca &  Barbuy ( 2017 )  with formation time-scale of 1 Gyr, for several radii, are o v erplotted: r < 0. 5 kpc (solid lines), 0. 5 < 1 kpc (dotted lines), 1 < r < 2 kpc (short-dashed lines), 2 < r < 3 kpc (long-dashed lines). stars older than 8 Gyr,  Ryde et al. ( 2016 ) and Siqueira-Mello et al. ( 2016 ). Fig. 10  shows that a typical star formation time-scale of 1 Gyr (the chemodynamical model with ν SF = 1 Gyr -1 ) also explains the Si and Ca abundances found in the bulge.
C O N C L U S I O N S We have selected 58 stars from the bulge sample of  Queiroz et al. ( 2021 )  based on APOGEE with characteristics to belong to the spheroidal pressure-supported bulge. For this sample, we have analysed lines of C, N, O, alpha-elements Mg, Si, Ca, and neutroncapture element Ce. Our fits to the Mg, Si, and Ca abundances from the original APOGEE results using the ASPCAP software appeared to be in good agreement. We recomputed abundances for C, N, O, and Ce, assuming the spectroscopic non-calibrated stellar parameters from APOGEE DR17. We report differences in abundances of these elements. We compare the abundances of these elements to literature data for bulge stars and chemodynamical models by Fria ca & Barbuy ( 2017 ) -see also  Barbuy et al. ( 2018a ). These comparisons show compatibility of the abundances of the sample stars with literature and models for Mg, Si, and Ca in which a pressure supported component (spheroidal bulge) formed on a very short time-scale (below 1 Gyr).
D. G. gratefully acknowledges support from the  ANID   BASAL  project  ACE210002. D. G. also acknowledges finan-cial support from the  Direcci ón de Investigaci ón y Desarrollo de la Universidad de La Serena  through the  Programa de Incentivo a la Investigaci ón de Acad émicos (PIA-DIDULS). The work of V. M. P. is supported by  NOIRLab , which is managed by the  Association of Universities for Research in Astronomy (AURA)  under a cooperative agreement with the  National Science Foundation. MZ was funded by  ANID FONDECYT Regular   1191505 ,  ANID Millennium Institute of Astrophysics (MAS)  under grant  ICN12 009 , the  ANID BASAL Center for Astrophysics and Associated Technologies (CATA)  through grants  AFB170002 ,  ACE210002  and  FB210003. DM gratefully acknowledges support by the  ANID BASAL  projects  ACE210002  and  FB210003  and by F ondec yt Project No.
This landscape is often regarded to be shaped by an underlying gene regulatory network (GRN), which appears as a powerful abstraction for describing interactions between genes through their proteins production. The construction of GRNs from literature being a very time-consuming and labor intensive process, and sometimes impossible due to the limitation of our current knowledge, their automated reconstruction from large datasets has become a classic task in systems biology  [44]. This task is notoriously difficult, in particular when dealing with singlecell transcriptomics, the bursty synthesis of mRNAs  [58] ,  [35]  giving rise to highly variable and non-Gaussian expression data  [26]. The methods that are used cover a wide range of statistical and modeling tools  [1] , including the analysis of stochastic models of gene expression  [16] ,  [5].
In the latter case, expression datasets are identified to independent samples of the time-varying distribution describing the process. GRN inference can then be seen as the reconstruction of the most-likely GRN from a set of partial observations of independent realizations of the model. To our knowledge, the use of metastability for performing such reverse engineering has not been studied yet. The main contribution of this article is to derive an efficient algorithm for linking the landscape analysis of a mechanistic model of gene expression using metastability, to the most-likely associated GRN parameters. In Section 1, we are going to present a mechanistic model of gene expression developed in  [16] , which describes single-cell dynamics associated to a GRN with transcriptional bursting.
With a methodology similar to the one developed in  [51] , we perform its reduction into a discrete coarse-grained model on a limited number of metastable basins. We deduce from this reduction an approximation of the time-dependent proteins distribution of the original model, presented in Section 2, which appears as a mixture of Gamma distributions. We develop in Section 3 a statistical method for linking the parameters of such mixture to the GRN parameters of the model. In Section 4, we extend this method for estimating GRN parameters from scRNA-seq timestamped datasets. We show in Section 5 the accuracy of the method for in silico datasets simulated from the mechanistic model with various networks of different sizes. Finally, we discuss more precisely in Section 6 the interpretation of the method in terms of landscape, its applicability to real datasets, and we highlight its similarity with a machine learning approach.
The parameter β i represents the basal activity of gene i, and each parameter θ ij encodes the interaction j → i. Every function k θ on,i is then comprised between two positive constants k 0,i < k 1,i , and ∂ P j k θ on,i has the sign of θ ij. This sigmoidal form for the bursts rate functions can be interpreted as a simplification of the mechanistic form used in  [16] ,  [5]. We point out that assuming that this model is able to reproduce real single-cell datasets, we make the underlying hypothesis, which is implicit in most GRN inference methods, that, contrary to its state, the GRN structure is not modified under the action of some hidden variables. One should note that our interaction model is an approximation of the underlying biochemical cascade reactions, and that the genes we modeled need not to be transcription factors.
Indeed, any GRN θ is associated for any value of ε to an unknown stationary distribution, that we denote ûθ ε , which converges to a sum of Dirac on the attractors when ε → 0  [36]. Then we see that the Gamma mixture ûε associated to the same ε converges to the same sum of Dirac when ε → 0, as the mean of each Gamma does not depend on ε while the variance is proportional to ε: we have |û θ ε -ûε | → 0 as ε → 0. When the noise ε is not negligible, the problem is more complex and finding a theoretical bound for the quantity |û θ ε -ûε |, depending on ε, is beyond the scope of this paper. Nevertheless, we see in Figure  3  that the Wasserstein distance between the empirical distribution associated to a proteins dataset simulated from the mechanistic model (3) and the Gamma mixture distribution  (6)  is much smaller than the distance between the same dataset and a mixture of normal distribution (fitted with a Gaussian Mixture Model).
Altogether the reduction and methods described above establish a formal basis for the definition of a simplified epigenetic landscape given a GRN, under the form of a mixture model. As the mixture parameters seem possible to obtain from a single-cell dataset (see Section 4. 2), it would be fruitful to use the same formalism to assess the inverse problem of inferring the most likely GRN, given an (experimentally-determined) cell distribution in the gene expression space, a notoriously difficult task  [39] ,  [16]. 3 Method for inferring a GRN from a mixture distribution In this section, we discuss the problem of identifying the most-likely GRN which is associated to an approximate landscape described by the Gamma mixture distribution of the form  (6). We identify such landscape to the set of parameters: α = µ(z), k z,i d 1,i , c i z∈Z,i=1,••• ,n.
Clustering From the statistical model  (10) , the first step consists in inferring for each timepoint t k the set of parameters α k M = {µ t k (z), α z , c} z∈Z , where c = k off k off +s 0 ∈ R n and for every z ∈ Z, α z = kz d 0 ∈ R n. The number of genes being potentially large, a minimization of the maximum likelihood function f on the gene expression matrix M t k ∈ M nt k ,n (N), defined by f (M t k ) = nt k j=1 z∈Z µ t k (z) n i=1 NB k z,i d 0,i , k off ,i k off ,i + s 0,i (x k j ) , even using an EM algorithm or an other variational method, would be very uncertain, such algorithms being known to be trapped into local minima  [10]. The number of components of the mixture, i. e the size of the set Z, remains also unknown. A good method for fitting α k M would be to apply a Bayesian procedure like a Reversible-Jump Monte-Carlo Markov Chain (RJMCMC) algorithm  [42]  to the whole dataset.
Finally, the regression step consists in solving successively, for k = 1, • • • , l, the problem θ k = arg min θ R k (θ, α k M ) + λ|θ -θ k-1 |, (16) where θ 0 is defined as the null matrix of size n, |θ -θ k-1 | is given by (  15 ), λ is a penalization coefficient and the function R is defined for every θ and α k M = (µ t k (z)) z∈Z , (α z ) z∈Z by: R k (θ, α k M ) = n i=1 z∈Z µ t k (z) σ θ i ( α z α 1 ) - α z,i α 1,i 2. ( 17 ) The procedure is illustrated in Figure  5. For the applications presented in Section 5, the value of the coefficient λ has been calibrated in order to be optimal for various datasets simulated from randomly-generated tree-like networks, like those studied in Section 5. 6. Interestingly, the cost function (  17 ) can be obtained as the particular case of a more general method for linking a GRN to a set of mixture parameters in the case where single-cell proteomic data were available, and that mRNAs are seen as a proxy for the proteins level.
We present this method in Appendix C. We underline the fact that the function σ θ , defined in (2), depends only on the GRN θ and αz α 1 , which can be directly estimated from scRNA-seq data, using the method described in Section 4. 2. 1 and the relations  (12). Thus, we do not need to make any assumption on the value of the hyperparameters of the model, not even the ratio d 0 d 1 as it was the case in  [15]. This is due to the fact that we do not infer the protein distribution but only the values of its main modes, which are completely characterized by the mean of the proteins distribution within each basin, and do not depend on its variance. The two-steps method presented in Sections 4. 2. 1 and 4. 2. 2, that we call CARDAMOM (Cell types Analysis from scRna-seq Data Achieved from a Mixture MOdel), shows good results for the simulated datasets that are presented in Section 5.
Back to the model: consistency of the algorithm and verifications The inference method presented above corresponds to the calibration of the mechanistic model (1) driven by burst rates functions of the form (2). Recalling that GRN inference is generally not an identifiable problem, many networks are likely to reproduce a dataset. Beyond the topology of the network that we expect to be well reconstructed by the algorithm, a GRN given by CARDAMOM should be considered as accurate if it allows to reproduce the datasets used for the inference, seen as partial observations of the time-varying distribution associated to the mechanistic model. It would then be natural to use the quantitative parameters inferred with CARDAMOM to simulate new snapshot data that we could compare to the original data.
Simulating data with stimulus The simulations of the model (1) are performed with the python package HARISSA presented in  [15] , which is based on an efficient thinning method for simulating the model  (1). For reproducing in vitro differentiation processes, we first simulate every network until reaching a first equilibrium before t = 0. At t = 0, we introduce a new gene called "stimulus", the proteins level of which is artificially maintained at the value of 1. This value corresponds to to its highest possible attractor: indeed, the proteins scaling  (11)  implies that no coordinates of the attractors can be greater than 1 in the gene expression space. Proteins being much less noisy than mRNA (i. e d 0,i d 1,i 1), the protein level associated to a gene is expected to be close to the attractor associated to the basin it belongs to, and forcing a gene to be fully expressed in the model  (1)  with the scaling (  11 ) is equivalent to force its value to 1.
Inference of a toggle-switch network from a static expression data sampled under the long time distribution We first evaluated the ability of CARDAMOM to find the parameters of a 4-genes network consisting in 2 independent toggle-switch of the same form than the one described in Table  1. To this aim, we simulate 10 dataset by sampling 50 independent cells at 2 timepoints 0 and 20h. In that case, the stimulus has no effect on the system. Inference is then independently performed for every datasets and the mean of all the significant values obtained for each network is presented in Figure  6 , compared to the network used for the simulations. We observe in Figure  6  that the algorithm allows to reconstruct very well the topology of the network, detecting which pairs of genes are in relation together and the sign of the interactions.
Inference of a 4-genes network with branching and feedback loop from timestamped data We now evaluate CARDAMOM on the 4-gene network described in Figure  5  and Table  2. Although such a small network may appear very simple, it already has some interesting features (branching, feedback loop with inhibition) and is interesting for inference. To this aim, we simulate 10 datasets by sampling independent cells at 10 time points t = 0, 2, 4, 6, 8, 11, 13, 15, 17, and 20h, with 50 cells per timepoint (then 500 cells in total for each dataset). Inference is independently performed for every dataset and the results are merged into receiver operating characteristic (ROC) and precision-recall (PR) curves that are shown in Figure  7. CARDAMOM turns out to reconstruct very efficiently the topology of the network.
We also compare in Figure  10  the two temporal dynamics of the two complete datasets, one used for the inference and one simulated with the inferred network. The data are projected altogether with the method umap  [27]  and we show on each subfigure the cells corresponding to one of the datasets. These cells are then classified depending on the time where the measures are realized after applying the stimulus. The results suggest that the model has been relatively well calibrated through the inference procedure.
We underline that this approximation does not mean that the data always exhibits strong bimodality. It only means that the fact that they are sometimes observed far from their extremum, which is expected at least during the period of transition after applying the stimulus, does not mean that these transitory states are stable. The simplification consisting in imposing a number of clusters equal to 2 for each gene would not affect the efficiency of the algorithm for the 4-genes network studied in Section 5. 3. Indeed, we observe in Figure  8 (A) that the RJMCMC algorithm applied to each gene of the network for every datasets simulated from this network usually computes a number of clusters equal to 2, and in Figure  8 (B) that imposing a number of cluster greater than 2 makes slightly decreases the efficiency of the algorithm.
Application to larger networks We now consider the cases of tree-like activation networks of 5, 10, 20, 50 and 100 genes. For each case, we simulate ten datasets corresponding to 10 random networks of the same size, that sampled from the uniform distribution over trees rooted in the stimulus  [15]. All datasets contain the same 10 timepoints t = 0, 2, 5, 8, 11, 13, 16, 19, 22, and 25h, with 100 cells per timepoint (then 1000 cells in total for each dataset). Inference is then independently performed for all dataset with CARDAMOM and the state-of-the-art method GENIE3  [19]. The results are presented in terms of Area Under Precision-Recall curves (AUPR), in Figure  11A ). CAR-DAMOM turns out to reconstruct more efficiently the topology of the network than GENIE3. The strong decrease in performance when the number of genes gets higher is due not only to the lack of data per timepoints in regards to the number of interactions, but also to the choice of the timepoints.
Indeed, a sequence of timepoints that is too coarse to catch the dynamics would lead to a lack of accuracy in the inference, and a sequence which is too tight would often be too short, leading to miss the activity of some genes. This appears to be the case from 20 genes in Figure  11A ). For every size of network, an average runtime is obtained after inferring the 10 datasets associated to the 10 tree-like networks, on a 16-GB RAM, 2,4 GHz Intel Core i5 computer. We see in Figure  11B ) that for both algorithms the computational speed increases linearly with respect to the number of genes, with a slope which is significantly lower for CARDAMOM. These results show that the algorithm is suitable for realistic number of genes and cells.
Discussion and prospects In this discussion, after clarifying in Section 6. 1 the importance of the clustering step of CAR-DAMOM, we go deeper in Section 6. 2 into the link between the mixture approximation and the popular notion of Waddington landscape. Then, we discuss in Section 6. 3 the applicability of our method to real single-cell datasets. Finally, we show in Section 6. 4 that the regression step of the method highlight the link existing between the mechanistic model and a neural network, and in what way this analogy paves the way for a more general method using adaptive bursts rates functions.
Importance of the clustering The attentive reader should have remarked that the function  (17)  which is optimized in the second step of CARDAMOM appears very simple and close to a generalized linear model, which is surprising for an analysis based on a mechanistic model. It is then natural to ask whether the clustering step remains important or if it is possible to reduce the method to the regression step by identifying each cell to a distinct attractor without losing an important accuracy (the clustering step would then simply consists in rescaling each count of a gene i by the parameters c i ). We compare in Figure  12  this coarse simplification with CARDAMOM on the datasets used in Section 5. 3. We observe that such simplification generates a significant decrease in accuracy, confirming that the method cannot be reduced to such generalized regression.
Applicability of CARDAMOM to real single-cell datasets The algorithm CARDAMOM has been shown to reconstruct accurately a most-likely GRN from timestamped in silico datasets, through the estimation of the metastable parameters associated to these data. It is then natural to ask whether the method could deal with real datasets or not. The method for estimating the metastable parameters should be well suitable for real datasets, provided that the Negative Binomial mixture approximation of mRNAs is accurate, and we have seen in Section 5. 6 that the algorithm is suitable for realistic number of genes and cells. The main difficulty that could appear would then be related to the decrease in the performance of the algorithm with the number of modes appearing in the sample. Indeed, we recall that the regression step aims to infer a GRN matrix of size n × n from the values of the functions k on,i at each attractor.
To go further, we recall that the mathematical analysis beyond CARDAMOM lays on the point of view that cell differentiation processes can be coarsely reduced into a discrete process on a limited number of cell types, i. e that the main ingredients for characterizing stochasticity are the frequency modes describing the cell types and the random transitions between them, which is commonly accepted since  [18]. Such transitions have been for example recently proposed as the basis for facilitating the concomitant maintenance of transcriptional plasticity and stem cell robustness  [55] : in this case, the authors had proposed a phenomenological view of the transition dynamics between states, and our work may typically connect this cellular plasticity to an underlying GRN dynamics. This connection between a GRN and associated cell states should also be used for the quantitative modelling of stochastic state transitions underlying the generation of diversity in cancer cells  [56] ,  [14] , or for investigating transitions between the large diversity of clusters that can be observed in human ES cell differentiation datasets  [29].