Globalizing comparative political science research on Southeast Asia

From my point of view, the scholar of Asia today is a 'global scholar', embedded not only in his/her area studies specialization, but also in constant intellectual dialogue with others in the same discipline, but focused on other Area studies fields. In practice, by dint of the complexity of the 'object', we are studying to becoming interdisciplinary. In short, the scholar of Asia today not only has greater savoir, s/he has an enhanced savoir faire: in a virtuous circle of scientific dialogue, greater interdisciplinarity leads to today's Asia scholar making a more substantial contribution to his/her 'home' discipline. Below, I shall argue the case for political science.

David Camroux

study of Southeast Asia in which he underlines the importance of state institutions and political parties in determining equitable outcomes, an argument applicable outside of the region. [START_REF] Kuhonta | The Institutional Imperative: The Politics of Equitable Development in Southeast Asia[END_REF] All three place emphasis on the role of external forces and particularly on the immediate post-World War II period of decolonization as being a watershed in the path dependencies they describe. This insight is also crucial in Ja Ian Chong's, a Singaporean based academic's, nuanced comparative study of state formation in Indonesia and Thailand, as well as China. [START_REF] Ian | External Intervention and the Politics of State Formation: China, Indonesia and Thailand, 1893-1952[END_REF] Three of the above volumes have been published by Cambridge University Press, which has emerged as the most exciting publisher in the area of Southeast Asian politics. CUP has also just published an edited volume by one of the three editors of the 2008 Stanford publication, Erik Martinez Kuhonta, who with Allen Hickens brings together a representative group of younger political scientists to discuss political parties in Asia in a conceptually rich way. [START_REF]Party System Institutionalization in Asia: Democracies, Autocracies, and the Shadows of the Past[END_REF] CUP has also published two textbooks by more senior scholars, Jacques Bertrand and Bruce Gilley, both of which attest to the mainstreaming of the new scholarship in Southeast Asian politics. [START_REF] Bertrand | Political Change in Southeast Asia[END_REF] This is also the case in two comparative overviews of the state of democracy throughout Asia, published by Palgrave Macmillan in 2010 and 2014. [START_REF] Paul | Obstacles to Democratization in Southeast Asia: A Study of the Nation State, Regional and Global Order[END_REF] Routledge, alongside Palgrave Macmillan, the main British non-university press publishers, have contributed to this flowering of political science research on Southeast Asia by bringing to a much wider audience the work of a number of dynamic European political scientists: Marco Bunte, Aurel Croissant, Dirk Tomsa and Andreas Ufen. [START_REF]The Crisis of Democratic Governance in Southeast Asia[END_REF] Their work, like that of a previously German-based academic Mark Thompson, 17 is significant both in its comparative dimension and also in its delving into European political science traditions. Routledge has also has been very active in publishing the work of Asian 18 and Australian scholars. [START_REF] Rich | Parties and Parliaments in Southeast Asia: Non-partisan chambers in Indonesia, the Philippines and Thailand[END_REF] In crossing disciplinary boundaries, anthropology has become even more central to the research agenda of political scientists as evidenced in an edited volume published in 2012 on Southeast Asian Perspectives on Power. [START_REF]Southeast Asian Perspectives on Power[END_REF] Terence Lee at the National University of Singapore draws on anthropology in his study of military responses to popular protests in Asia, and in doing so provides theoretical insights that are salient beyond the Asian cases he examines. [START_REF] Lee | Defect or Defend: Military Responses to Popular Protests in Authoritarian Asia[END_REF] A mere seven years following the Stanford volume, to simply provide a measure of how far we have come, allow me to mention the fourth, very recently published volume in the Routledge Handbook collection devoted to Southeast Asia. [START_REF]Routledge. The other three handbooks deal with history, politics and economics[END_REF] Edited by William Case, a prolific writer from the intermediate generation, [START_REF] Case | Politics in Southeast Asia: Democracy or Less, London: Routledge Curzon, represented a pioneering effort in comparative analysis of Southeast Asian polities. Director of the Southeast Asia Research Centre at the City University of Hong Kong[END_REF] and entitled the Routledge Handbook of Democratization in Southeast Asia, it is a representative display of 'the state of the art', theoretically rigorous and resolutely comparative.

To conclude, in order to provide for the busy reader a very recent, readily available, example that supports the central argument of this brief overview: an extended article by Marcus Mietzner on the presidency of Joko Widido (Jokowi) in Indonesia published online by the East West Center. [START_REF] Mietzner | Reinventing Asian Populism: Jokowi's Rise, Democracy, and Political Contestation in Indonesia[END_REF] Mietzner, a Canberra-based European scholar who along with Edward Aspinall and Greg Fealy has made the ANU a mecca for the study of Indonesian politics. In this article Mietzner, by postulating a concept of 'technocratic populism' to describe Jokowi's praxis, and juxtaposing this with other theoretical arguments drawn essentially from Latin America, demonstrates how the vibrant Southeast Asian political experience requires us to rethink a number of assumptions and interpretations based on observations in other regions of the world. IN LESS THAN A DECADE two interrelated developments have impacted on both the discipline of political science and the field of Southeast Asian Studies. On the one hand a new generation of scholars in North America, Europe, Australasia and Southeast Asia itself has emerged to revitalize research in political science in relation to the ten ASEAN member countries. On the other, the Southeast Asian experience has come to generate theory rather than, as was previously largely the case, being an object on which theory and theoretical concepts from outside were brought to bear.

David

Of course this latter development had previously occurred in other fields. Consider the work of two intellectual giants in their respective disciplines: Benedict Anderson in history and James Scott in anthropology. Anderson's broad path-breaking work on Indonesia led directly to his seminal study on nationalism globally, Imagined Communities; first published in 1983, one that continues to feed into all debates on this subject. [START_REF] Anderson | Imagined Communities: Reflections on the Origin and Spread of Nationalism[END_REF] Similarly, James Scott's Weapons of the Weak, based on a study of peasant farmers in Malaysia, led to a general rethinking of power relations in rural societies globally. [START_REF] Scott | Weapons of the Weak: Everyday Forms of Peasant Resistance[END_REF] While working on themes for his marvelous study of communities in the uplands of continental Southeast Asia, Living Without a State, Scott published an inspiring trans-national study on the State and on State capacity. [START_REF] Scott | The Art of Not Being Governed: An Anarchist History of Upland Southeast Asia[END_REF] In the recent past one can cite only two authors in political science/political economy who have had an impact on theorization in a global context. Jeffrey Winter's important study of oligarchy in Indonesia, Power in Motion, [START_REF] Winters | Power in Motion: Capital Mobility and the Indonesian State[END_REF] published in 1996, led later to his widely praised comparative study of varieties of oligarchy, in which, as a heuristic devise, he makes an international comparison between warring, ruling, sultanistic and civil oligarchies. [START_REF] Winters | Oligarchy[END_REF] Muthiah Alagappa, in his conceptual introductory chapters to the five volumes he edited, sought to draw from the broader Asian experience to break new ground in thinking on subjects such as political legitimacy, civil-military relations, the study of international relations and processes of democratization. [START_REF]Civil Society and Political Change in Asia: Expanding and Contracting Democratic Space[END_REF] These five studies were published in the 'Contemporary Issues in Asia and the Pacific' series at Stanford University Press, a cooperative undertaking with the East West Center, with Alagappa himself the series editor. It therefore comes as no surprise that, with John Sidel and Geoffrey White succeeding him as editors, that my first 'bookend' title featuring the new generation of political science scholars, was published within this series.

The publication in 2008 of an edited volume, Southeast Asia in Political Science, demonstrated that comparative political science research on Southeast Asia had come of age. [START_REF]Southeast Asia in Political Science: Theory, Region, and Qualitative Analysis[END_REF] The eleven chapters by a new generation of scholars rapidly making their mark (as well as one by Don Emmerson, a pioneer in the field) are all theoretically rich and avowedly comparative. 8 The three editors of this volume -Dan Slater, Erik Martinez Kohunta and Tuong Vu -have since continued to make a significant contribution to political science research. Dan Slater's innovative study of state formation in Southeast Asia has undoubtedly added substantially to the general comparative theoretical literature on State formation. [START_REF] Slater | Ordering Power: Contentious Politics and Leviathans in Southeast Asia[END_REF] His typology of institutional outcomes and their causes has relevance beyond the Southeast Asian case studies he conducted, as does the link between contentious and conciliatory politics and authoritarian and democratic outcomes. Like his co-editor Tuong Vu, Slater emphasizes the importance of critical historical junctures in creating political path dependencies. Tuong Vu's broad canvas study comparing emerging, or non-emerging, development trajectories in two Southeast Asian countries, Indonesia and Vietnam, with those in two Northeast Asian nations, South Korea and China, is theoretically sophisticated. [START_REF] Vu | Paths to Development: South Korea[END_REF] He posits a link between patterns of elite alignment and elite-mass engagement and the emergence, or nonemergence of developmental states. Martinez Kuhonta returns to the subject of development in his own broad

  Camroux is Honorary Associate Professor and Senior Researcher at Sciences Po (Centre for International Studies), Paris. (david.camroux@sciencespo.fr).


Collaboration over Competition: How US Institutions Should Adapt to the Growth of the Asian Higher Ed Market evolllution.com/opinions/collaboration-competition-institutions-adapt-growth-asian-higher-ed-market

Asian universities are growing in strength and stature, and it would behoove leading US institutions to collaborate with these growing institutions rather than ignore their rise.

There is no doubt about Asia's presence in the international education landscape. No region has undergone such a profound transformation of its higher education in the past decades, the dynamics being closely interconnected with the unprecedented economic growth. From Northeast to South, most Asian countries have recently experienced major social and demographic changes and institutional reforms, often leading to greater political stability. The growth of a large middle class, coupled with the openness driven by economic imperatives, has pushed national governments to massively expand access to higher education while bringing local institutions to the level set by international standards.

What makes the specificity of the evolution of the higher education sector in Asia in the past 15 years is the central role that national governments are asking universities to play in the broader economic development. Because economic growth is, today more than ever, related to knowledge production and advanced skills, Asian governments and citizens expect higher education institutions to create the conditions for the development of their countries and to train the future generation to be innovative and creative and to pursue sustainable growth. Hans de Witt and Philip Altbach call the phenomenon "nationalization of internationalization," a clear shift that is not unique to Asia, though the means and pace of reforms are.

Across a variety of political and institutional settings, historical legacies, and local constraints, the reform of higher education in Asia has been largely state-driven and selectively publicly funded. On the contrary, in other regions, the state gradually withdrew from investments in higher education, and encouraged private investments and market competition while not solving inclusiveness and equal opportunities concerns. What all higher education reforms in Asia have in common is that they took inspiration from overseas, while mainly addressing local issues. China, Japan, South Korea and Taiwan, as well as ASEAN members and India have regularly referenced Western models and practices-most often those of the Unites States, but also Europe. Gaining visibility and becoming more compatible with the rest of the world was crucial to the effectiveness of all the local policies aiming at encouraging foreign investments, reversing the brain drain, attracting international students, creating interregional student mobility and transferring know-how through international joint initiatives. International rankings have captured some of the positive effects of these reforms that are allowing the most established national universities to emerge internationally and acquire global visibility and reputation. According to different ranking methods (THE, QS and ARWU), Asia has this year between 17 and 36 universities among the world's top 200. Yet, the transformation is much deeper than brand development effects: The gross enrollment ratio (GER) increased in all the countries, in order to serve the aspirations of a growing and wealthier domestic population and to boost local talent pools. Data collected by the World Bank reveal that not only countries like Japan and Korea have some of the highest GER in the world (61 percent and 98 percent), but also that the rest of Asia reduced the gap in a spectacular way in less than 15 years. China moved from 5 percent to 27 percent, India from 17 percent to 25 percent, Indonesia from 12 percent to 32 percent, Vietnam from 3 percent to 25 percent, Thailand from 32 percent to 51 percent, Malaysia from 11 percent to 37 percent, just to mention some of them.

These numbers are the result of proactive public policies. Funding was key, of course, even though public expenditure in tertiary education as a percentage of the GDP per capita did not necessarily grow, but in most cases remained stable or decreased. Almost all Asian governments, with the exception of India, have adopted an elitist approach, concentrating strong financial support to a small number of universities, whose potential to match the international standard was high. In addition, like in many Western countries, reforms brought financial decentralization, corporatization, autonomy and competitive allocation of resources. They often unlocked bureaucratic restrictions, encouraged partnerships and introduced some elements of market competitions, even in economies that are largely state driven.

The impressive trajectory of Asian countries is not undermining the leading role of U.S. universities, particularly the research-based schools, which have become the reference of a world-class university. Even if the new metrics introduced by international rankings since 2003 (i.e. employer reputation, research achievements of the faculty, and quality of services) give a comparative advantage to universities able to invest in faculty and brand development, assuring consistent quality remains a key issue in most Asian countries, especially with the growth in the number of private institutions that are not publicly regulated.

What we observe though the lens of Asia is an evolution of higher education similar to other business sectors, particularly in the service industry: Newcomers are challenging the incumbents, as they are often better positioned strategically to adapt to the new needs of the students, being more flexible and faster in implementing innovation. The past 15 years show that Asia has become a laboratory with great ambitions, with remarkable capacity to change, and with the joint support of national governments and the economic major players, both interested in innovation and skills.

The trend is reinforced by a growing volatility of students' international trajectories. Even if the US and the UK are still the top destinations for international students, flows now grow in a number of other directions, including from the west towards Asia. Students aspire at accumulating multiple experiences, including studying in several countries. Asia has strong assets, linked to its economic growth and the employment and business opportunities Studying in Asia is no longer limited to those who are passionate about oriental language and civilization. Some U.S. universities have understood this very early and created branch campuses or platforms in Asia, students today have direct access to information. They can choose and organize their studying experience in Asia without the mediation of a foreign university. The best and the brightest do not need to go though the traditional study abroad. Programs like the Schwartzman scholar program or the NUS Yale College of liberal arts are appealing directly to them.

Parochialism is today the major threat to U.S. universities, whose established research excellence prevents them from noticing the changes in the whole playing field. Working in a more collaborative manner with Asian universities via dual degrees, joint teaching and combined curricula might be the only way not to miss the boat.


Thibault A Etienne

Charlotte Roux

Eugenio Cinquemani

Laurence Girbal

Muriel Cocaign-Bousquet

Delphine Ropers
email: delphine.ropers@inria.fr

A nonlinear mixed-effects approach for the mechanistic interpretation of time-series transcriptomics data

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Temporal omics data provide unprecedented views of cellular inner workings. They account for the molecular responses of living organisms to perturbations, resulting from regulatory mechanisms that ensure adaptation and survival. Statistical modelling tools are commonly used to analyse and interpret these data, from the clustering of gene expression profiles according to their similarity or characteristics to network inference (for reviews, [START_REF] Ding | Analysis of time-series regulatory networks[END_REF][START_REF] Grigorov | Analysis of time course omics datasets[END_REF][START_REF] Lu | High-dimensional ODEs coupled with mixed-effects modeling techniques for dynamic gene regulatory network identification[END_REF][START_REF] Noor | Biological insights through omics data integration[END_REF][START_REF] Straube | A linear mixed model spline framework for analysing time course 'omics' data[END_REF]. Achieving mechanistic insight requires moving from data-driven to mechanistic modelling approaches, by including prior knowledge of the relations between molecular responses and cell components [START_REF] Androulakis | Analysis of time-series gene expression data: Methods, challenges, and opportunities[END_REF][START_REF] Noor | Biological insights through omics data integration[END_REF]. Mechanistic model parameters have a biological interpretation that provides insight on molecular mechanisms at stake. Mechanistic modelling is, however, a daunting task and examples of these models are poorly abundant in the literature [START_REF] Noor | Biological insights through omics data integration[END_REF]. Main issues arise at the time of estimating the parameters due to model and data complexity. Combinations of technological and practical limitations result in small numbers of biological and technical replicates that can be used and of time points that can be measured. Appropriate data normalisation is required to disentangle sources of noise, resulting from measurement errors and biological process heterogeneity, within and across studies [START_REF] Liu | Normalization methods for the analysis of unbalanced transcriptome data: a review[END_REF]. The consequence of these shortcomings has been a coarse-graining of the models to help with estimation, at the price of an incomplete interpretation of the data: much biological insight remains buried.

The aim of this article is to develop a method for the estimation of large-scale mechanistic model parameters from dynamical transcriptomics data, addressing the above issues. We extend the nonlinear mixed-effects (NLME) modelling framework that allows to study variability between individuals in a population. Frequently used in pharmokinetics-pharmocodynamics studies of the variability of patient responses to drug absorption and elimination, NLME models have been successfully applied in other scientific fields as well (epidemiology, agriculture...) and, more recently, in cell biology to study cell-to-cell variability [START_REF] Almquist | A nonlinear mixed effects approach for modeling the cell-to-cell variability of Mig1 dynamics in yeast[END_REF][START_REF] Davidian | Nonlinear models for repeated measurement data: an overview and update[END_REF][START_REF] Lavielle | Mixed effects models for the population approach: models, tasks, methods and tools[END_REF][START_REF] Lindstrom | Nonlinear mixed effects models for repeated measures data[END_REF][START_REF] Llamosi | What Population Reveals about Individual Cell Identity: Single-Cell Parameter Estimation of Models of Gene Expression in Yeast[END_REF][START_REF] Loos | A hierarchical, data-driven approach to modeling single-cell populations predicts latent causes of cell-to-cell variability[END_REF][START_REF] Oddi | A nonlinear mixed-effects modeling approach for ecological data: Using temporal dynamics of vegetation moisture as an example[END_REF]. In NLME modelling, the same individual response model is used to describe a physical or biological phenomenon observed for each individual and a statistical model describes why the model parameters vary between individuals. The information shared between individuals results in better estimates for both the individual and population parameters [START_REF] Almquist | A nonlinear mixed effects approach for modeling the cell-to-cell variability of Mig1 dynamics in yeast[END_REF][START_REF] González-Vargas | Validation methods for population models of gene expression dynamics[END_REF].

Unlike traditional NLME approaches, we model "mRNA-to-mRNA" variability, in an effort to identify biochemical principles behind the biological variability of mRNA concentrations. The idea is that the same cellular machineries determine mRNA intracellular levels, whose mechanism of action can be described by the same individual response model for each mRNA. For instance, in bacteria, a unique RNA polymerase transcribes all mRNAs and a single multiprotein complex, degradosme, is responsible for bulk mRNA degradation. Messenger RNAs can thus be regarded as forming a population of individuals subjected to the same cellular machineries. The variability between individual mRNAs reflects the intervention of additional regulatory mechanisms adjusting transcription or degradation to environmental conditions. In this study, we will focus on mRNA concentrations averaged over populations of cells, because many transcriptomics data are obtained in these conditions. We will show that NLME modelling is a suitable framework to analyse large, dynamical microarray and RNA-Seq datasets: 1) it separates the technical and biological variability of high-throughput data; 2) it leverages the limited number of replicates and time points by imposing a parameter distribution over the population; (3) it allows a knowledge-based interpretation of the datasets by incorporating mechanistic details; (4) it empowers a posteriori analyses of the biological variations between mRNAs to identify regulatory targets. Our method is a versatile approach, generalisable to other case studies.

We illustrate our approach by applying it to the study of bacterial mRNA decay, a key cellular process to adjust gene expression to environment [START_REF] Dressaire | Role of mRNA stability during bacterial adaptation[END_REF][START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF]for review, Carpousis 2007). Degradation of bacterial mRNAs by the degradosome leads to short-lived messengers, typically between 3 and 8 minutes in Escherichia coli, up to 50 to 100 minutes in some extreme cases [START_REF] Bernstein | Global analysis of mRNA decay and abundance in Escherichia coli at single-gene resolution using two-color fluorescent DNA microarrays[END_REF][START_REF] Chen | Genome-wide study of mRNA degradation and transcript elongation in Escherichia coli[END_REF][START_REF] Dressaire | PNPase is involved in the coordination of mRNA degradation and expression in stationary phase cells of Escherichia coli[END_REF][START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF]. The half-lives are generally determined following transcriptional arrest, by means of transcriptomics experiments monitoring at the genome-wide level the residual mRNA concentrations along time. The degradation profile of each mRNA reflects the activity of degradosome and regulatory mechanisms operational at the time of drug addition. It is classically fitted by simple exponential models, assuming a first-order reaction rate and independent degradation reactions, from which is then determined mRNA half life [START_REF] Bernstein | Global analysis of mRNA decay and abundance in Escherichia coli at single-gene resolution using two-color fluorescent DNA microarrays[END_REF][START_REF] Chen | Genome-wide study of mRNA degradation and transcript elongation in Escherichia coli[END_REF][START_REF] Laguerre | Large-scale measurement of mRNA degradation in Escherichia coli: To delay or not to delay[END_REF][START_REF] Moffitt | Spatial organization shapes the turnover of a bacterial transcriptome[END_REF]. In general, the degradation profiles are relatively similar, stemming from the fact that mRNAs are degraded in bulk by the same machinery, but differ at the level of mRNA abundance prior to drug addition and in their curvature. To mechanistically characterise this variability, we adapted a genome-wide model of mRNA degradation described in [START_REF] Etienne | Competitive effects in bacterial mRNA decay[END_REF]. It explicitly relates the degradation profiles with the degradation machinery and additional regulatory mechanisms involving, for instance, small RNAs and ribonucleoproteins [START_REF] Etienne | Competitive effects in bacterial mRNA decay[END_REF]. Using our NLME framework, we estimated mechanistic model parameters from a microarray and a RNA-Seq dataset on Escherichia coli bacteria growing in exponential phase [START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF]Roux et al., in preparation). The method uncovers mRNAs known to be post-transcriptionally regulated, notably by small RNAs, such as mRNAs coding for enzymes in the central carbon metabolism and for ribosomal proteins.

Modelling approach

Fig. 1 gives an overview of our method, which takes as input preprocessed microarray or RNA-Seq time-series data. A mechanistic model relates their time-dependent evolution with cellular components and regulatory mechanisms, but model parameters are unknown. We model "mRNA-to-mRNA" variability using the NLME framework to estimate population and individual mRNA parameters. Biological variation between individual parameter values can be interpreted in relation to regulatory mechanisms.

Transcriptomics data

The two transcriptomics datasets used in this study were obtained in the same conditions [START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF]Roux et al., in preparation). They quantify the residual mRNA concentration over time of E. coli bacteria grown exponentially at 0.63 h -1 in minimal medium supplemented with glucose until antibiotic treatment by the transcriptional inhibitor rifampicin. Three samples were taken right before the addition of rifampicin, and nine after. Measurements for each of the 4254 mRNAs at each time point consist of eight microarray probe Figure 1: Nonlinear-mixed effects approach for the mechanistic interpretation of time-series transcriptomics data. The framework uses preprocessed microarray or RNA-Seq data, obtained from samples collected at time point t j with j = {0, ..., J} for cellular genes i = {1, ..., I}. They quantify the abundance mi (t j ) of each mRNA. Probe intensities r = {1, ..., R} are available for microarray data. An optional smoothing step enables to remove their variability (Section 3.2). NLME modelling using a mechanistic model gives estimates of population-and individual-mRNA parameters. The variability of the later has biological interpretation.

intensities corrected for the background or raw read counts for the RNA-Seq (Supplementary Section S1). Lowly expressed genes with signal intensities or read counts close to the background were filtered out. This leaves us with the same 2809 mRNAs in each dataset.

Mechanistic modelling

A variety of models can be used in NLME modelling. Here we use an ordinary-differential equation (ODE) model to describe mRNA concentrations averaged over a population of cells. We have shown in [START_REF] Etienne | Competitive effects in bacterial mRNA decay[END_REF] that bacterial mRNA degradation can be assimilated to a macro-reaction catalysed by the degradosome, which we assume to follow Michaelis-Menten kinetics to reproduce the monotonous decrease of mRNA concentrations following transcription arrest. Owing to the fact that microarray and RNA-Seq experiments quantify mRNA abundances, we reformulate the model to describe the time-dependent evolution of mRNA quantities (see Supplementary Section S2 for details):

dmi(t) dt = - Vm × mi(t) Kmi + mi(t) V , with: mi(t0) = mi(τi), ( 1
)
where m i (t) is the abundance of mRNA i (expressed in AU or count) at time t (min) and V is the cell volume (in milligram of cell dry weight, mgDW). Vm, the maximal velocity of the reaction (in UA or count/mgDW/min), is identical for all mRNAs. Km i is the Michaelis constant of a given mRNA i (in UA or count/mgDW). Addition of rifampicin at time 0 inhibits the initiation of transcription. Degradation starts after a short delay τ i (in min), during which RNA polymerase already engaged in the elongation of transcription completes mRNA synthesis. Km i is an inverse measure of the affinity of degradosome for its target mRNA. For instance, mRNAs with small Km i value are degraded faster than mRNAs with an average Km. This value depends on the sequence and structure characteristics of the mRNA, but also on regulatory factors such as small RNAs or RNA-binding proteins, which can either facilitate or block the binding of degradosome to mRNAs.

Formulation of the NLME problem

The fact that the same machinery degrades the bulk of mRNAs allows to consider a population approach based on NLME modelling for the estimation problem. We consider that the bulk mRNAs form a joint population, wherein the degradation kinetics of each individual mRNA can be described by the same response model, described in Eq. 1. Model parameters for each individual mRNA are drawn from the parameter distribution of the population. Similarly to NLME approaches applied to single-cell data [START_REF] Dharmarajan | A simple and flexible computational framework for inferring sources of heterogeneity from single-cell dynamics[END_REF][START_REF] Llamosi | What Population Reveals about Individual Cell Identity: Single-Cell Parameter Estimation of Models of Gene Expression in Yeast[END_REF], the objective is to harness differences between mRNA levels in order to improve the estimation of population distribution parameters and of individual parameters. Concretely, we describe the time-series data with the following measurement model:

mi,t j = f (mi(tj), Φi) + g(f (mi(tj), Φi), θ) × ϵi,j. (2)
Here, mi,tj is the measured concentration of mRNA i = 1, . . . , N at time t j , ϵ i,j ∼ N (0, 1) is the residual error, and the function g the residual error model with a vector of noise parameters θ. Function f is the solution to the ODE model in Eq. 1. The vector of parameters Φ i = {m i (t 0 ), τ i , Vm i , Km i } is drawn from a probabilistic distribution h characterised bya vector Φ pop corresponding to the typical value of model parameters, also called fixed effects, and a variance-covariance matrix of the random effects Ω describing the variability between individual mRNAs. Here we will consider that h is a lognormal distribution to constrain individual parameters to positive values:

Φi ∼ LogN (Φpop, Ω) ⇔ log(Φi) = log(Φpop) + ηi , ( 3
)
with η i ∼ N (0, Ω). Additional levels of variability can be considered, for instance in Section 3.2, to explicit technical variability between microarray probes [START_REF] Lavielle | An improved SAEM algorithm for maximum likelihood estimation in mixtures of non linear mixed effects models[END_REF][START_REF] Panhard | Extension of the SAEM algorithm for nonlinear mixed models with 2 levels of random effects[END_REF]. This requires to introduce a new vector of parameters, Φ i,r , for individual mRNA i and probe r:

Φi,r ∼ LogN (Φpop,r, Ω) ⇔ log(Φi,r) = log(Φpop,r) + ηi + ηi,r (4)
with η i ∼ N (0, Σ), η i,r ∼ N (0, Γ), and Ω = Σ + Γ [START_REF] Panhard | Extension of the SAEM algorithm for nonlinear mixed models with 2 levels of random effects[END_REF]. Σ is the variancecovariance matrix for inter-individual variability (between mRNAs) and Γ the variance-covariance matrix for intra-individual variability (between probes of a given mRNA).

The residual error model g in Eq. 2 should be adapted to the transcriptomics data. Measurement variance is considered to be proportional to probe intensities or read counts [START_REF] Meacham | Identification and correction of systematic error in high-throughput sequence data[END_REF][START_REF] Weng | Rosetta error model for gene expression analysis[END_REF]. Systematic errors in transcriptomics data include non-zero background intensity level and hybridisation bias

Vmpop [AU (count)•mgDW -1 •min -1 ] 1.77e +3 ± 10.50 3.64e +3 ± 25 1.15e +5 ± 834 Kmpop [AU (count)•mgDW -1 ] 4.28e +3 ± 75.1 7.01e +3 ± 108 1.7e +5 ± 2.19e +3 Catalytic efficiency [min -1 ] 0.41 0.52 0.68 m(to)pop [AU (count)]
10.9 ± 0.24 8.57 ± 0.20 426 ± 13.9

τpop [min] 0.33 ± 1.66e -3 0.52 ± 1.7e -3 0.74 ± 2.1e -3 ωKm [AU (count)•mgDW -1 ] 0.80 ± 1.37e -3 σKm = 0.71 ± 1.19e -3 0.50 ± 9.63e -4 γKm = 0.16 ± 1.88e -4 ω m(to) [AU (count)] 1.17 ± 1.58e -3 σ m(to) = 1.25 ± 1.66e -3 1.72 ± 2.33e -3 γ m(to) = 0.55 ± 2.87e -4 ωτ [min] 1.65 ± 4.90e -3 στ = 1.50 ± 3.1e -3 1.22 ± 2.46e -3 γτ = 0.14 ± 6.3e -4 a [AU (count)] 0.84 ± 5e -4 0.24 ± 8.75e -5 3.82 ± 7.26e -3 b [-] 0.06 ± 5.61e -5 0.15 ± 3.79e -5 0.24 ± 1.55e -4
Condition number 15.0 10.6 5.9

related to the binding efficiency of microarray probes for instance, or bias related to the library size, transcript length, and GC content in the case of RNA-Seq data. Data preprocessing allows to remove part of the systematic error, through background correction of probe intensities and normalisation to reduce the variability between measurements and time points in the case of dynamical data. We use a combined error model to take into account the technical errors of RNA-Seq and microarrays:

g(f (mi(tj), Φi), θ) = f (mi(tj), Φi)) × b + a) (5)
where θ = {a, b}. The multiplicative error parameter b is typically related to the mean of probe binding efficiencies and a is the basal measurement error [START_REF] Meacham | Identification and correction of systematic error in high-throughput sequence data[END_REF][START_REF] Weng | Rosetta error model for gene expression analysis[END_REF]. These error parameters will be estimated together with the population parameters (Φ pop , Σ, Γ, and Ω), from which, together with the data, are then derived the specific parameters for individual mRNAs.

Estimation of mechanistic model parameters from microarray data

Estimation from smoothed data

Can we estimate mechanistic model parameters from a single dynamical microarray dataset using NLME modelling? In this section, we address the question, by working on the most favourable scenario. We preprocess the background-corrected microarray probe intensities to take out part of the systematic measurement error. Among the possible approaches [START_REF] Laguerre | Large-scale measurement of mRNA degradation in Escherichia coli: To delay or not to delay[END_REF][START_REF] Oh | Temporal dynamic methods for bulk RNA-Seq time series data[END_REF][START_REF] Ritchie | limma powers differential expression analyses for RNA-sequencing and microarray studies[END_REF], we have chosen spline smoothing to adjust for variability between the eight probes and the twelve time points for each mRNA i [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF][START_REF] Ritchie | limma powers differential expression analyses for RNA-sequencing and microarray studies[END_REF] (Supplementary Section S1.1). We obtain one degradation profile per mRNA and a total of 2809 smoothed degradation profiles to be used for estimation of the fixed and random effects of model parameters. The later are defined as follows: mRNAs being degraded at the same maximal rate, Vm is set to a fixed value, while Km, the delay before degradation τ , and initial concentrations at time t 0 , m(t 0 ), are allowed to vary between mRNAs. This gives:

       log(Vmi) = Vmpop log(τi) = log(τpop) + η τ i log(m(t0)i) = log(m(t0)pop) + η m(t 0 ) i log(Kmi) = log(Kmpop) + η Km i with η ϕ i ∼ N (0, σ 2 ϕ ) for ϕ = {τ, m(t 0 ), Km}.
Estimation of NLME model parameters was run with the Monolix software in 4997 CPU seconds (see Supplementary Section S3 for details). As shown in Table 1 (second column), the population parameters are estimated with high certainty, with a small standard error for most of them and small error model parameters. The good estimation result reflects the size of the dataset (2809 degradation profiles) used for estimating parameter distributions. Our mechanistic model is sufficiently complex to explain the data, according to the normality of the error residuals with a mean of 0 (Supplementary Section S4.1). The Visual Predicted Check in Supplementary Fig. S5 shows that model simulations with the estimated parameters are able to reproduce the central trend and variability of the degradation profiles. This is illustrated in Fig. 2(a) for mRNA satP. The predicted degradation profile using the estimated individual parameters correctly fits the data. This mRNA has a lower initial concentration at time t 0 , and a slower degradation kinetics than the typical mRNA, representative of the population, whose degradation profile is shown as a black line.

Michaelis-Menten equations are known to pose practical identifiability problems, notably to obtain unique estimates of Vm and Km, while their ratio is identifiable (e.g., [START_REF] Choi | Beyond the Michaelis-Menten equation: Accurate and efficient estimation of enzyme kinetic parameters[END_REF][START_REF] Holmberg | On the practical identifiability of microbial growth models incorporating michaelis-menten type nonlinearities[END_REF][START_REF] Stroberg | On the estimation errors of KM and V from time-course experiments using the michaelis-menten equation[END_REF]. In the NLME framework, the fact that individual parameters are tied by the population distribution can render otherwise intractable identification problems feasible [START_REF] Lavielle | What do we mean by identifiability in mixed effects models?[END_REF]. Our population parameters seem to be practically identifiable, since they have low standard errors and the estimation algorithm converges to almost the same solution when started several times from random initial conditions (Supplementary Fig. S8). In addition, the condition number corresponding to the ratio of eigenvalues of the variance-covariance matrix of the estimates is satisfactory (Table 1). Correlations between parameters obtained from this matrix are shown in Fig. 2(d). All of them are statistically significant (p-value < 2.2• × 10 -6 with a t-test). Most likely, our parametrisation choice contributed to the small correlation of 0.35 between population parameters Vm pop and Km pop (with individuals set to a fixed Vm value, while Km i varies), which can be seen in the multistart estimation results in Supplementary Fig. S8. The two parameters, not just their ratio, are informative. We observe correlations between other estimates, though. For instance, between m(t 0 ) pop and Km pop (-0.66) and between Ω τ and τ (-0.61). It may reflect residual identifiability issues (e.g., a bimodal distribution of τ pop would have been more appropriate, see Supplementary Fig. S11 for η τ i ), as well as a biological phenomenon: mRNA abundance is a known determinant of mRNA stability [START_REF] Esquerré | Genome-wide investigation of mRNA lifetime determinants in Escherichia coli cells cultured at different growth rates[END_REF][START_REF] Nouaille | The stability of an mRNA is influenced by its concentration: a potential physical mechanism to regulate gene expression[END_REF], which in our model is largely determined by the Km parameter. This may explain the correlation between m(t 0 ) pop and Km pop . In the same line, mRNAs in E. coli are mainly regulated at the transcriptional level [START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF]. Abundant mRNAs are thus actively transcribed and on-going transcription at the time of rifampicin addition is likely to delay the onset of degradation, hence the small correlation of 0.33 between abundance and delay at the population level.

Are the estimated parameters biologically relevant? We cannot say for Vm pop and Km pop whose nonstandard units prevent their comparison to literature data. Their ratio in standard units, 0.41 min -1 in Table 1, represents the catalytic efficiency of degradosome and mRNA turnover. It indicates that the typical individual mRNA is degraded every 2.4 min, which agrees well with the median half life of 2.8 min determined for the same dataset in [START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF]. The τ pop value in Table 1 corresponds to an average delay of about 20s before the onset of degradation, in agreement with the expected delay of rifampicin action <30s observed experimentally (e.g. [START_REF] Chen | Genome-wide study of mRNA degradation and transcript elongation in Escherichia coli[END_REF][START_REF] Pato | Residual rna synthesis in Escherichia coli after inhibition of initiation of transcription by rifampicin[END_REF]. E. coli has many long mRNAs such as polycistronic mRNAs, which often delay the onset of degradation [START_REF] Chen | Genome-wide study of mRNA degradation and transcript elongation in Escherichia coli[END_REF]. This most likely explains the relatively large value of the random effect for parameter τ (ω τ = 1.65 min).

The largest population parameter for random effect is ω m(t0) , which means that differences in mRNA concentrations explain a large part of the variability of the degradation profiles. Km values vary as well, but to a lesser extent. They reflect variations of the affinity of the degradosome for its target mRNAs and possible post-transcriptional regulations. This will be further studied in Section 5. We conclude from these results that NLME modelling is a suitable framework for mechanistic parameter estimation from dynamical transcriptomics data, at least when data is preprocessed. Despite the fact that only one dataset is used for estimation, the numerous individuals in the population share sufficient information to obtain good inference results.

Estimation from microarray probe intensities

In this section, we complexify the estimation problem. Since mRNA abundance is quantified by means of eight probes varying in their binding efficiency and mRNA complementarity, we have a population of 22,472 probes and as many degradation profiles, which could be useful information to improve parameter inference. However, probe intensities are background corrected but not further preprocessed. Can the NLME model handle richer, but noisier and larger data, at the price of an increased computational cost? Technical variability between probes can be seen for mRNA satP in Fig. 2(b) (black points). While the curvature of the eight degradation profiles is similar, probe intensities are variable at time t 0 , all the more so due to the use of three measurements for each probe at this time point.

We introduce a new random effect in our NLME model, η i,r , to describe technical variability between probes specific to a given mRNA i: with η ϕ i ∼ N (0, σ 2 ϕ ) and η ϕ i,r ∼ N (0, σ 2 ϕ ) for ϕ = {τ, m(t 0 ), Km}. As expected, the estimation problem is computationally more expensive, with a CPU time of 2.4•10 6 s. Accurate and robust estimation results were obtained with the microarray probe intensities, as shown in Supplementary Information S4 for the residual distributions, VPC, and the convergence of the multistart algorithm to the same solution. The condition number of the variance-covariance matrix of the estimates is smaller when probe intensities are used (Table 1) and correlations between population parameters are markedly reduced too (Fig. 2(e).

       log(Vmi) = Vmpop log(τ i) = log(τpop) + η τ i + η τ i,r log(m(t0)i) = log(m(t0)pop) + η m(t 0 ) i + η m(t 0 ) i,r log(Kmi) = log(Kmpop) + η Km i + η Km i,r
The γ parameter values in Table 1 (third column), indicative of the technical variability between probe intensities, remains smaller than the variability between mRNAs, which is biological variability. The largest random effect between probes and mRNAs, observed for γ m(t0) , reflects the important noise of probe intensities at time t0. Some of them have clearly a weaker binding efficiency, leading to an estimated initial concentration m(t 0 ) below that obtained with smoothed data. Interestingly, the standard deviation of random effects between mRNAs, σ Km , σ m(t0) , and σ τ , which corresponds to biological variability, is similar in both cases. They lead to consistent random effects between mRNAs for both datasets (Supplementary Fig. S11). For instance, mRNAs with weak affinities have high Km random effects in both datasets, and conversely. In the case of the random effects for initial concentrations, strong consistency between estimates from smoothed and nonsmoothed data is achieved when considering the two levels of variability (Supplementary Fig. S12), which indicates that information is gained when explicitly accounting for fluctuations between probes in the model.

The only changes, concerning error model parameters, are indicative of data preprocessing being performed within the NLME framework. The absence of data smoothing is reflected in the doubling of parameter b value, associated with binding efficiency bias. To the contrary, explicit consideration of the technical variability has decreased the basal measurement error represented by parameter a. The parameter value might represent residual systematic bias related to acquisition and quantification of array images for instance, or a systematic hybridisation bias. We note a small divergence between the Vm pop and Km pop values obtained from smoothed and nonsmoothed data. We cannot exclude that estimates from smoothed data were biased by the preprocessing step performed outside of the NLME framework in Section 3.1.

Overall, these results suggest that mechanistic parameter estimation from noisy and large microarray data is possible with the NLME framework. It allows accurate statistical inference from (almost) raw data through a good separation of the technical and biological variability.

Estimation of mechanistic model parameters from RNA-Seq data

RNA sequencing has become a popular approach for the quantification of gene expression. It provides a higher coverage and greater resolution than microarrays and former DNA sequencing approaches, and alleviate hybridisation bias [START_REF] Marguerat | RNA-seq: from technology to biology[END_REF][START_REF] Hrdlickova | RNA-Seq methods for transcriptome analysis[END_REF][START_REF] Todd | The power and promise of rna-seq in ecology and evolution[END_REF]. Thus, with the wealth of RNA-Seq data being generated, it is interesting to inquire whether they are also prone to mechanistic interpretation. The difficulty is that we have only one noisy time series per mRNA for parameter estimation, as illustrated in Fig. 2(c) for mRNA satP.

NLME parametrisation and error modelling are the same as described in Section 3.1. There are no technical replicates, thus technical random effect has been omitted, and we consider only one level of variability, between mRNAs. Population parameters are estimated accurately, with a low standard error (Table 1, fourth column). The mechanistic model is able to explain count data and describe their variability, and the estimation results are robust, as shown in Supplementary Section S4. RNA-Seq and microarray units are different, which limits the comparison of estimation results. The Vm pop and Km pop values lead to a somewhat higher, but similar mRNA turnover, in the range of those measured experimentally (Table 1) [START_REF] Bernstein | Global analysis of mRNA decay and abundance in Escherichia coli at single-gene resolution using two-color fluorescent DNA microarrays[END_REF][START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF]. The value of the average delay, τ pop , remains close to the experimental observations [START_REF] Chen | Genome-wide study of mRNA degradation and transcript elongation in Escherichia coli[END_REF][START_REF] Pato | Residual rna synthesis in Escherichia coli after inhibition of initiation of transcription by rifampicin[END_REF].

In the absence of technical random effect, part of the technical variability is accounted for by the error model parameters, which are higher (Table 1). While we have opted for the same combined error model for microarray and RNA-Seq data for the sake of comparison, the correlation between a and b values indicates that a simplified model would be more appropriate (-0.53, Fig. 2(f)). Apart from Vm pop and Km pop that are more correlated in this estimation (0.52), other correlations are weak and the condition number of the variance-covariance matrix is the smallest of the three estimations (Table 1). Hence, even in the absence of technical replicate, the estimated biological random effects are consistent with that obtained from another type of data (Supplementary Fig. S11).

Altogether, inference of mechanistic model parameters from RNA-Seq data is also possible with the NLME framework, despite data noise and the absence of technical replicate. Mechanistic models provide physiologicallyinterpretable parameters, which will be analysed below. In the context of NLME models, in particular, we would like to mechanistically interpret the random effects, because they represent the effect of underlying regulatory mechanisms.

Identification of regulatory mechanisms of mRNA degradation

Classically, transcriptomics data are used in differential expression analyses, conducted on the fold change of expression levels characterising gene expression variability between conditions. Our NLME framework gives us access to random effects that represent individual variations from the mean, mode, or median. These can be interpreted as fold changes and can be used directly in Gene Set Enrichment Analysis (GSEA) [START_REF] Yu | clusterProfiler: an R package for comparing biological themes among gene clusters[END_REF].

Here we will perform GSEA on one specific model parameter, the Michaelis constant. The random effect associated with this parameter, η Km i , subsums fluctuations of mRNA sequence and structure characteristics affecting the binding of the degradation machinery, as well as the effect of regulatory factors facilitating or inhibiting this binding (see Supplementary Eq. S2). Activatory effects decrease the value of individual parameter Km i and thus, the corresponding random effect η Km i can take a negative value, while inhibitions increasing the Km lead to higher, positive random effects. Non-regulated mRNAs are expected to have small random effects, close to zero. They are already well described by population parameters.

Fig. 3 shows the results of a GSEA conducted on the parameters estimated from the smoothed microarray dataset. It reproduces the η Km i distributions of mRNAs in the GSEA enriched Gene Ontology (GO) categories. One category, organic acid processes related to carbon metabolism, is strongly post-transcriptionally regulated. The category includes genes involved in glycolysis and the Krebs cycle, whose half life is known to be controlled at high growth rate [START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF]see Roux et al. 2021 for review). It is interesting to note that mRNAs involved in the organic acid catabolic process category tend to be stabilised. This agrees with the stabilisation of mRNAs involved in the COG category Carbohydrate transport and metabolism in [START_REF] Esquerré | Genome-wide investigation of mRNA lifetime determinants in Escherichia coli cells cultured at different growth rates[END_REF]. It may reflect the need to accumulate more of the corresponding enzymes to sustain high growth rate on glucose, which is the experimental condition studied. On the contrary, the negative random effect of Km for mRNAs in organic acid biosynthetic processes indicates that they are less needed in these growth conditions. These processes are involved in fast adaptation to fluctuating environments. Alleviating this post-transcriptional control can trigger a quicker cellular response through rapid accumulation of the corresponding enzymes, as mRNA instability and the relatively low cost of transcription compared to protein synthesis are determinant for fast changes of gene expression [START_REF] Pérez-Ortín | Homeostasis in the central dogma of molecular biology: the importance of mRNA instability[END_REF]. The five other categories shown in Fig. 3, encompassing aminoacid biosynthesis, cell motility, and translation, are known to include numerous post-transcriptionally regulated genes [START_REF] Modi | Functional characterization of bacterial sRNAs using a network biology approach[END_REF][START_REF] Wagner | Small RNAs in bacteria and archaea: who they are, what they do, and how they do it[END_REF]. According to our study, the degradation of these genes is stimulated, although it can be compensated for by high transcription rates. For instance, genes in the for the smoothed microarray data (Supplementary Section S6). Distributions are coloured according to their adjusted p-value translation and ribosome assembly categories, notably genes coding for ribosomal proteins, are known to be actively transcribed at high growth rate, in particular in these experimental conditions [START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF]. In our study, this is reflected by the high random effect obtained for the initial concentrations of these mRNAs.

The above results show that our method is able to retrieve known post-transcriptionally regulated genes. We can now make hypotheses on the regulatory mechanisms involved. Half of the genes in each category shown in Fig. 3 are known and/or predicted targets of small RNAs [START_REF] Modi | Functional characterization of bacterial sRNAs using a network biology approach[END_REF][START_REF] Wright | CopraRNA and IntaRNA: predicting small RNA targets, networks and interaction domains[END_REF]. The main ones are cited on the right side of the figure. They facilitate or block degradosome binding to the target mRNA, acting alone or in complex with a RNA-binding protein such as Hfq or ProQ (for reviews, [START_REF] Holmqvist | RNA-binding activity and regulatory functions of the emerging sRNA-binding protein ProQ[END_REF][START_REF] Modi | Functional characterization of bacterial sRNAs using a network biology approach[END_REF][START_REF] Reyer | Kinetic modeling reveals additional regulation at co-transcriptional level by post-transcriptional sRNA regulators[END_REF][START_REF] Santiago-Frangos | Hfq chaperone brings speed dating to bacterial sRNA[END_REF]. For instance, the small RNA gcvB is the master sRNA regulator of amino-acid metabolism and transport in a wide range of Gram-negative bacteria [START_REF] Urbanowski | The gcvB gene encodes a small untranslated RNA involved in expression of the dipeptide and oligopeptide transport systems in Escherichia coli[END_REF]Miyakoshi et al. 2022 and references therein). It represses more than 50 genes in E. coli at the post-transcriptional level [START_REF] Miyakoshi | Mining RNA-seq data reveals the massive regulon of GcvB small RNA and its physiological significance in maintaining amino acid homeostasis in Escherichia coli[END_REF], in agreement with the negative random effect of individual Km parameters observed for the category cellular amino acid biosynthetic process in Fig. 3.

Small RNAs have emerged as major post-transcriptional regulators that affect translation initiation and/or mRNA stability, both negatively and positively [START_REF] Wagner | Small RNAs in bacteria and archaea: who they are, what they do, and how they do it[END_REF]. Yet their targets have not be all identified. For instance in Fig. 3, sole half of the genes in the different GO biological process categories have been identified as post-transcriptionally regulated. Our approach has thus allowed to identify new target genes.

Discussion and conclusion

The approach developed in this paper responds to the pressing need of mechanistically interpreting highthroughput biological data. To that aim, we have extended the application of NLME models to analyse "mRNA-to-mRNA" variability in a population of cells, by considering that mRNAs form a joint population submitted to the same cellular machinery, degradosome in our study. This allows to infer multidimensional parameter distributions describing the population from dynamical transcriptomics data, and then to derive specific parameters for individual mRNAs. We have shown that the framework allows consistent interpretation of heterogeneous datasets with different levels of technical noise. It factors out technical variability and incorporates mechanistic details, thus allowing to leverage biological variability for the identification of regulatory mechanisms and their targets.

Such mechanistic interpretation of dynamical omics data is still rare in the literature, because statistical inference of nonlinear models is challenging. One notable exception is the estimation of a model of cancerrelated signalling pathways from exome-and transcriptome-sequencing data [START_REF] Fröhlich | Efficient parameter estimation enables the prediction of drug response using a mechanistic pan-cancer pathway model[END_REF]. The 4100 unknown parameters of this large ODE model have been accurately estimated from data comprising more than 6,900 experimental conditions. Imposing a parameter distribution over the population in our study leverages the limited number of experimental observations despite a larger number of unknown parameters. This is a remarkable result, as this implies that reliable estimation from less comprehensive high-throughput datasets is possible.

This study provides the first example of the same NLME framework applied to different types of data, and to such a large amount of longitudinal data. We have shown that the framework is flexible and can be used to analyse data of different types, with different levels of technical noise, from raw RNA-Seq data to nonsmoothed and smoothed microarray probe intensities. The framework is flexible and can be extended to other case studies, with different data, models, conditions, organisms, in single cells or populations of cells... provided adaptation of the model for the source of variability and of the residual error model. In the context of mRNA degradation, the possibility to analyse a large number of data is a clear improvement over classical studies, wherein determination of reliable mRNA half lives is often restricted to about 20% of cellular mRNAs (e.g., [START_REF] Chen | Genome-wide study of mRNA degradation and transcript elongation in Escherichia coli[END_REF]. Useful information is lost, which makes it difficult to draw general conclusions from a subset of data.

There is an advantage to integrate additional data, for instance to improve parameter identifiability. In our study, the addition of literature information such as the value of the maximal velocity could solve the small discrepancies observed for Vm pop and Km pop (Section 3). However, integration of heterogeneous datasets is easier if they are expressed in compatible units. This would require an absolute quantification of transcriptome to obtain data in standard international units.

Current limitation to the integration of larger datasets is essentially computational. Multiplying by eight the dimension of the data and adding one level of variability in Section 3.2 has increased by ∼ 50 the CPU time. The problem arises when the dependency structure of the individual parameters is complex. Efforts to overcome this limitation have started, notably with the development of f-SAEM, a fast mixing sampler for the maximum likelihood estimation of NLME population parameters [START_REF] Karimi | f-SAEM: A fast stochastic approximation of the EM algorithm for nonlinear mixed effects models[END_REF], or of the platform Pumas with parallelization capabilities for simulation and estimation [START_REF] Rackauckas | Accelerated predictive healthcare analytics with Pumas, a high performance pharmaceutical modeling and simulation platform[END_REF]. Another means of addressing this issue is data reduction, for instance through nonparametric smoothing of aggregated individuals with similar gene expression profiles [START_REF] Lu | High-dimensional ODEs coupled with mixed-effects modeling techniques for dynamic gene regulatory network identification[END_REF].

Data preprocessing

Transcriptomics data used in this study come from microarray and RNA-Seq experiments. Their preprocessing was carried out using the R software (R Core Team, 2021).

Microarray data

The dataset is published in [START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF]. It corresponds to E. coli K12 MG1655 cells grown exponentially on glucose at 0.63h -1 prior to addition of the transcriptional inhibitor rifampicin. Cultures were sampled before and after drug addition and total RNA was extracted. Two µg of total RNA were converted to cDNA and subsequently labelled and hybridised to E. coli gene expression arrays. Two probesets of eight different probes were used per gene. Averaging the signal between probesets provided us with the intensity of eight different probes per mRNA. Probe intensities were then corrected for background using random probesets [START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF]. Rather than applying an intrareplicate quantile normalisation to the intensities of the eight probes used per mRNA, as done in [START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF], we corrected the variation between probes and between time points through spline smoothing [START_REF] Wang | Shape-restricted regression splines with R package splines2[END_REF]. This allowed to obtain a single time series per gene. An example of spline smoothing is provided for mRNA thrC in Fig. 4.

An additional filtration step was performed to retain only 2809 genes sufficiently expressed to allow for the observation of the degradation kinetics. Probe intensity for these genes is at least twice the background at time 0.

Abundance

RNA-Seq data

The strain and the experimental conditions are the same as in Section 6.1.1. Total RNA (1 µg) was ribodepleted using the riboPOOLs kit and used to construct a sequencing library using the Ion Total RNA-Seq Kit v2 and the Ion Xpress RNA-Seq Barcode Kit. Reads were mapped with the tmap package onto the E. coli genome (version NC_000913.3, GenBank) and counted with HTSeq-count ("intersection non empty" mode) (Roux et al., in preparation). We use directly the raw data of the 2809 genes selected in the microarray preprocessing step.

Mechanistic model of mRNA degradation

Simplification of mechanistic model of mRNA degradation

In [START_REF] Etienne | Competitive effects in bacterial mRNA decay[END_REF], the following mechanistic model describes the time-dependent evolution of the concentration c i (t) of each mRNA i, averaged over a population of cells subjected to treatment by rifampicin:

dci(t) dt = - kcat × E0 × ci(t) Km app i + E0 + ci(t) , ( 6
)
where kcat is the catalytic constant in min -1 , E0 the concentration of degradosome in µM, Km app i the apparent Michaelis constant of mRNA i in µM, and ci(t) the concentration of mRNA i in the same units. Km app i is said apparent because it reflects the effect of activators or inhibitors on the kinetics of the degradation reaction. For instance:

Km app i = Km t i 1 +
x small RN A K small RN A + K regulatory protein x regulatory protein + ... , (7) where the x's denote the concentration of regulatory factors and the K's, their dissociation constant. Km t i is the "true" Michaelis constant of mRNA i, which solely depends on the mRNA features. In this example, the small RNA has an inhibitory effect that decreases the affinity of the degradosome for its target mRNA, while the regulatory protein stimulates degradation.

For the purpose of this paper, we do not need to describe the enzyme concentration explicitly. We introduce new parameters: Vm denotes the maximal velocity of the reaction Vm = kcat×E0 (unit: µM•min -1 ), while Kmi = Km app i +E0 (unit: µM). This formulation leads to the classical Michaelis-Menten expression usually obtained after application of the standard quasi-steady-state assumption:

dci(t) dt = - Vm × ci(t) Kmi + ci(t) . ( 8
)
This is our starting model, which will adapt below so as to integrate the experimental data (Section S2).

Adaptation of mechanistic model to data

Mechanistic model with quantities While model variables and parameters are expressed in the International System (SI) of units, this is not the case of the microarray and RNA-Seq data. In addition, the latter two correspond to abundances of mRNAs.

The integration of data into the model can be performed in two different manners: either we transform the data to obtain concentrations, possibly in SI units if experimental design has allowed the absolute quantification of mRNA abundances, or we adapt the model to describe the time-dependent evolution of quantities in their units. We have chosen the latter approach in this study. This allows to estimate model parameters and initial concentrations easily, without the requirement of additional data to allow for data conversion. The inconvenient is that the estimates will not be expressed in SI units, which impedes their comparison with literature data. However, the objective of the study is to determine how the variability of individual parameters account for the variability of degradation profiles. We do not need SI units to reach such conclusions.

Let mi(t) denotes the quantity of mRNA i at time t. It is expressed in Arbitrary Units (AU) (microarrays), or in read counts (RNA-Seq). mi(t) is proportional to the cellular abundance ni(t) of mRNA i at time t averaged over the population of cells: mi(t) = α ni(t). ni(t) is expressed in µmol, while α is a conversion factor expressed in AU/µmol (microarrays) or in count/µmol (RNA-Seq). The concentration ci(t) is proportional to the quantity of mRNA i:

ci(t) = mi(t) V , (9) = α ni(t) V , ( 10
)
where V is the cell volume expressed in Liter (L). It is constant as rifampicin addition stops cell growth.

If we know α or if we measure the absolute abundance of mRNA i, ni(t), we can reformulate Eq. 8 as a function of ni(t). As this is not the case here, we replace ci(t) by its expression (9) in Eq. 8:

dmi(t) dt = - Vm × mi(t) Kmi + m i (t) V . ( 11
)
Explicit consideration of cell volume Analysing the degradation kinetics in various conditions may involve a change of volume between conditions, which implies to take into account the volume explicitly, as done in Eq. 11. Such information is not readily available in most cases, but alternative data can be used in place. In [START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF], for instance, the multiplication of the microarray probe intensities in AU, measured from a starting material of 2 µg of total RNAs, by the total RNA extraction yield (in µg total RNA per mg of dry cell weight) provides mRNA concentration values in UA/mgDW. The extraction yield Ey is inversely proportional to cell volume and thus can be used as proxy for the cell volume. It has been measured in the experimental conditions studied here: Ey = 45.2 µg/mgDW in E. coli cells growing exponentially at 0.63 h -1 [START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF], which gives V = 0.022 mgDW. We will use this value in our study. This changes the units of parameters Vm and Kmi to AU (or counts)/mgDW/min and AU (or counts)/mgDW, respectively.

Parameter estimation procedure

Initialisation of parameter values

Our data in AU or counts does not permit to use available literature values for initialising the unknown parameters. We used instead the mean of the degradation constants in [START_REF] Esquerré | Dual role of transcription and transcript stability in the regulation of gene expression in Escherichia coli cells cultured on glucose at different growth rates[END_REF], which is similar to the mRNA turnover or catalytic efficiency of degradation given by the ratio Vm/Km ≈ 0.5, from which initial values were given for the two parameters: Vm = 102 AU or count/gDW/min and Km = 180 AU or count/gDW based on several preliminary estimation tests. Standard deviations were initialised to 1. Parameters of the residual error model, a and b, were initialised to 1 and 0.3, respectively.

Parameter estimation with NLME modelling

Kinetic parameters were estimated using the SAEM-MCMC algorithm implemented in Monolix (Monolix version 2020R1. Antony, France: Lixoft SAS, 2020) using default settings indicated in the software documentation for auto-criteria of convergence at https://monolix.lixoft.com/ (maximum/minimun of the iterations, autocriteria, etc).

Estimation of population parameters

The estimation of the population parameters is the key task in nonlinear mixed effect (NLME) modelling. It was performed using the Stochastic Approximation Expectation-Maximization (SAEM) algorithm of Monolix. It is achieved in three steps:

1. The burn-in phase: it is the SAEM initialisation, wherein individual parameters are sampled from their conditional distribution by MCMC sampling using initial values of the population parameters. The number of iterations of this phase was set to five iterations.

2. The Exploratory phase, during which the parameter value at iteration k is built using information collected at that iteration only (it does not take into account the value of the parameter at the previous iteration). We set 500 iterations for the RNA-Seq and smoothed microarray data, and 1000 iterations for the nonsmoothed microarray data.

3. The smoothing phase. During this step the parameter value at iteration k is built using information collected at the previous iteration. We fixed 200 iterations for both datasets.

Estimation of individual parameters

The estimation algorithm for population parameters gives a rough estimate of the individual parameters already. However, it can be estimated by two more precise estimators: the conditional mode and the conditional distribution. We used them both.

Estimation of standard errors

The standard errors represent the uncertainty of the estimated population parameters. They were determined in Monolix through the calculation of the Fisher Information Matrix using stochastic approximation.

Assessing the quality of estimation

Distribution of residuals

Different types of residuals can be determined to detect potential misspecifications in the mechanistic model or the residual error model: the PWRES (Population Weighted RESiduals), the IWRES (Individual Weighted RESiduals), and the NPDEs (Normalised Prediction Distribution Errors). The distributions of the residuals are shown in Fig. 5 for the smoothed microarray data, Fig. 6 for the nonsmoothed microarray data, and Fig. 7 for the RNA-Seq data.

The following two tests indicate that the model correctly describe the experimental data, as residuals are symmetrically scattered around 0 and normally distributed (Shapiro-Wilk test).

Symmetry test:

Postprocessing of the estimation results

Functional enrichment analysis was performed with the R package clusterProfiler [START_REF] Yu | clusterProfiler: an R package for comparing biological themes among gene clusters[END_REF], by using the function gseGO. We chose the "ALL" ontology (biological process + molecular function + cellular component) and a pvalueCutoff equal to 0.05. We corrected the p-values with the Benjamini-Hochberg procedure. Functional groups possibly include mRNAs regulated at the post-transcriptional level by small RNAs. They can be identified by the prediction tool CopraRNA [START_REF] Wright | CopraRNA and IntaRNA: predicting small RNA targets, networks and interaction domains[END_REF].

Data availability

The Monolix script for estimation and output files are freely available at https://gitlab.inria.fr/tetienne/eccb_ script, together with the microarray dataset. The RNA-Seq dataset is being prepared for publication (Roux et al., in preparation) and will be made available on demand upon acceptance of the article.

Output files in .txt format include estimated population and individual parameters for each of the three estimations, together with the statistical results on model assessment (see the accompanying README file for details).

Figure 2 :

 2 Figure 2: Estimation of mechanistic model parameters from time-series transcriptomics data. Top panels: examples of fit for mRNA satP obtained with (a) smoothed microarray data, (b) microarray probe intensities, and (c) normalised RNA-Seq data. Continuous lines represent the predicted degradation profiles generated with estimated population parameters (black) and with individual parameters for mRNAs (red) and probe intensities (cyan). Bottom panels: correlation plots of the estimated parameters for (d) smoothed microarray data, (e) microarray probe intensities, (f) RNA-Seq data. All correlations are statistically significant (p-value < 2.2• × 10 -6 with a t-test). They are given according to the colour code shown below each plot.

  motility

Figure 3 :

 3 Figure 3: Gene set enrichment analysis conducted on random effects η Kmi

Figure 4 :

 4 Figure 4: Monitoring of thrC mRNA degradation by microarrays. The black dots represent the backgroundcorrected intensities of the eight probes at twelve time points, three time points before drug addition and nine after. The black line represents B-spline smoothing of the 8 probe intensities.

Figure 5 :

 5 Figure 5: Residual distributions for smoothed microarray data

Figure 7 :

 7 Figure 7: Residual distributions for RNA-Seq data

Figure 8 :

 8 Figure 8: Visual predictive check for the smoothed microarray data

Figure 10 :

 10 Figure 10: Visual predictive check for the RNA-Seq data

Figure 14 :

 14 Figure 14: Comparison of the random effects between mRNAs for τ (top panels), m(t 0 ) (middle), and Km (bottom), between smoothed and nonsmoothed microarray data (left column) and between nonsmoothed microarray and RNA-Seq data

Table 1 :

 1 Estimated mechanistic model parameters from smoothed microarray data, microarray probe intensities, and RNA-Seq data. The parameters of the distribution (mode and the logarithm of the standard deviation) are given for each NLME model parameter, together with their standard error. The catalytic efficiency is the ratio of Vm pop and Km pop . a and b are residual error model parameters. The condition number is the ratio of the eigenvalues of the variance-covariance matrix

	Data type and units:

Smoothed microarray data [AU] Microarray probe intensities

[AU]

Raw RNA-Seq data

[count]

Acknowledgements

The authors thank Inria MICROCOSME team members for numerous exchanges on this work. Aline Marguet (MICROCOSME) and Adeline Leclercq Samson (Univ. Grenoble Alpes) are more specifically thanked for discussions on mixed-effect modelling.

Funding

This work was supported by the French National Research Agency under project RIB-ECO [ANR-18-CE43-0010] and MEMIP [ANR-16-CE33-0018], as well as by the Institut National de Recherche pour l'Agriculture, l'alimentation et l'environnement (INRAE) and Inria, from which T.A.E. received a doctoral fellowship.

The Monolix script for estimation and output files are freely available at https://gitlab.inria. fr/tetienne/eccb_script, together with the microarray data. The RNA-Seq dataset is being prepared

Correlation matrix of the estimates

The inverse of the Fisher Approximation Matrix can be used to detect possible practical identifiability issues. It provides an approximation of the variance-covariance matrix of the estimated population parameters. A large condition number (the ratio of the largest and smallest eigenvalues of the matrix) is indicative of strong correlations between estimates. The (small) value of the condition number in our case is shown in Table 1 of the main text for the three estimations, while correlations are graphically represented in Fig. 2 of the main text.

Convergence assessment

Convergence of the estimation algorithm is another approach that can be used to detect identifiability issues [START_REF] Lavielle | What do we mean by identifiability in mixed effects models?[END_REF]. We started the algorithm from ten different initial conditions. Results are shown below in Figs. 111213. The algorithm converges to (almost) the same solution for each of the three estimations.