domain: sde math stat

text:
Chantal De Fouquet
email: chantal.de_fouquet@minesparis.psl.eu

Mathieu Le Coz

Xavier Freulon

Lea

Chantal De

Chantal De Fouquet

Léa Pannecoucke

Making kriging

Keywords: Groundwater modelling, Geostatistics, Radionuclide migration, Contamination

Introduction

Estimating the hydraulic head or concentrations of components in groundwater is necessary for drawing a map, locating "at best" new measurements by minimalizing redundancy [START_REF] Júnez-Ferreira | A geostatistical methodology for the optimal design of space-time hydraulic head monitoring networks and its application to the Valle de Querétaro aquifer[END_REF], or solving the inverse problem (Gomez-Hernandez et al., 1997;[START_REF] Angulo | Series Expansion Approach to the Inverse Problem[END_REF][START_REF] De Marsily | Modeling and Field Investigation in hydrogeology : a special volume in honor of Shlomo P. Neuman's 60th Birthday[END_REF]. However the head or concentration data are generally few and the kriging estimate of these variables with a usual variogram model is not necessarily consistent with the flow physics. Various methods have been proposed to constrain estimates by flow physics, a synthesis and non-exhaustive review of which is presented here.

As the head decreases in the flow direction, it is better to carry out the kriging using a non-stationary model [START_REF] Kumar | Optimal contour mapping of groundwater levels using universal kriging-a case study[END_REF]). An anisotropic variogram allows one to particularize the direction of flow at the macroscopic scale. Taking covariates into account allows for improvement of the modelling. Thus, [START_REF] Desbarats | On the kriging of water table elevations using collateral information from a digital elevation model[END_REF] introduce the topography as an external drift of the head, choosing the appropriate support of regularization of the digital elevation model. [START_REF] Varouchakis | Improving Kriging of Groundwater Level Data Using Non-Linear Normalizing Transformations -A field application[END_REF] make a non-linear estimation of the head (using a Gaussian anamorphosis or a modified Box-Cox transformation), with the necessary correction for the back transformation of the estimated values. On a real case, the cross-validation results are improved compared to the usual kriging.

Curvilinear coordinate transformation allows for bringing back the "natural coordinates" of flow, i.e. hydraulic head and stream functions [START_REF] Roth | Contribution de la géostatistique à la résolution du problème inverse en hydrogéologie[END_REF][START_REF] Rivest | Sparse data integration for the interpolation of concentration measurements using kriging in natural coordinates[END_REF]. In order to explicitly introduce flow physics, [START_REF] Brochu | A Simple Approach to Account for Radial Flow and Boundary Conditions When Kriging Hydraulic Head Fields for Confined Aquifers[END_REF] use the formalism of dual kriging to impose to the estimate the form of the analytical solution of a linear or radial flow calculated for a homogeneous medium. In [START_REF] Rivest | Hydraulic head field estimation using kriging with an external drift: A way to consider conceptual model information[END_REF] the covariate used as external drift is a « smooth » output of a finite element flow model calculated on an uniform transmissivity field, for which boundary conditions must therefore be specified. Tests on synthetic cases show that the accuracy of the estimate and its realism are improved. The method was applied to the real case of a dam.

On a synthetic case in steady state, [START_REF] Yang | Physics-Informed Kriging: A Physics-Informed Gaussian Process Regression Method for Data-Model Convergence[END_REF][START_REF] Yang | Physics-Informed Cokriging: A Gaussian-Process-Regression-Based Multifidelity Method for Data-Model Convergence[END_REF] reconstruct a field of concentrations by «physics-informed kriging» or «physics-informed cokriging». The method is based on a regression between observations which are considered as "high fidelity" data, and the average of outputs of flow and transport simulations which are considered as "low fidelity" data. The difference between kriging and cokriging lies in the covariance of the discrepancy between low fidelity data and observations. To improve the modelling, [START_REF] Yang | Physics-Informed Kriging: A Physics-Informed Gaussian Process Regression Method for Data-Model Convergence[END_REF][START_REF] Yang | Physics-Informed Cokriging: A Gaussian-Process-Regression-Based Multifidelity Method for Data-Model Convergence[END_REF] suggest a non-stationary model for this discrepancy.

In another approach, covariance models are constructed so as to be approximately or accurately consistent with flow equations. In the case of steady-state linear macroscopic flow and in the absence of recharge, the linearization of the differential equation linking perturbations of head and transmissivity provides an approximated solution. [START_REF] Hoeksema | An application of the geostatistical approach to the inverse problem in two-dimensional groundwater modelling[END_REF] supposed the head to be stationary. [START_REF] Dong | Kriging variables that satisfy the partial differential equation ΔZ=Y[END_REF][START_REF] Dong | Estimation géostatistique des phénomènes régis par des équations aux dérivées partielles (Geostatistical estimation of phenomena governed by partial differential equations[END_REF] gave a general solution within the framework of the intrinsic random function of order-k (IRF-k) theory and applied it to the inverse problem, by cokriging the transmissivity from transmissivity and head data. [START_REF] Hernandez | Conditioning mean steady state flow on hydraulic head and conductivity through geostatistical inversion[END_REF] introduced a source term in the equations for performing geostatistical inversion based on the pilot points method (de [START_REF] De Marsily | De l'identification des systèmes hydrogéologiques[END_REF][START_REF] De Marsily | Interpretation of interference tests in a well field using geostatistical techniques to fit the permeability distribution in a reservoir model[END_REF]. Zhang and Neuman (1996) generalized the analytical form of the simple and cross-covariance to velocity by crossing two flows with different main directions, for comparing longitudinal and transverse dispersion.

However the hypothesis of a uniform head gradient is often too strict. [START_REF] Roth | Combining geostatistics and flow simulators to identify transmissivity[END_REF] construct a variogram 𝛾(𝑥, 𝑥 ′) consistent with a steady state flow model with fixed aquifer geometry and boundary conditions. A set of geostatistical non conditional simulations of the conductivity are used as input of a flow simulation. The variogram 𝛾(𝑥, 𝑥 ′ ) is calculated, up to the ½ factor, as the mean square deviation of the head at pairs of points 𝑥 and 𝑥', computed on a set of flow simulation outputs. The numerical version of simple and cross-variograms (or covariances) between head and transmissivity is then used for cokriging these two variables. These experimental non-stationary variograms are used by Schwelle and Cirpka (2010) to estimate the head and the concentration by conditional expectation, calculated by the median (rather than the mean) on a set of simulations constrained by data from different variables. However, the conditional simulations of the transmissivity (or permeability) they use are obtained by an iterative inverse approach. [START_REF] Pannecoucke | Combining geostatistics and simulations of flow and transport to characterize contamination within the unsaturated zone[END_REF] generalize the approach to the vadose zone. Non conditional simulations of saturated hydraulic conductivity are used as input to flow-and-transport simulations of a non reactive contaminant. Boundary conditions can vary between the flow simulations, if these conditions are not precisely known. Applications to a wide variety of problems are presented in [START_REF] Pannecoucke | Combinaison de la géostatistique et des simulations à base physiqueapplication à la caractérisation de panaches de contaminants[END_REF] in two dimensions (2D) or 3D, for spatial or spatio-temporal estimation: delimitation of a contaminated zone at fixed contamination threshold and identification of the associated "uncertainty zone", backtracing of the unknown location of a contaminant source, comparison of the estimation with the results of a (real) 3D tracing experiment, etc. This paper presents original results of kriging with numerical covariance, following the approach of [START_REF] Pannecoucke | Combinaison de la géostatistique et des simulations à base physiqueapplication à la caractérisation de panaches de contaminants[END_REF], in order to take into account the heterogeneity of a contamination plume. However an anamorphosis (variable transformation) step is introduced, so that the contamination estimates are positive. This non-stationary anamorphosis is calculated on the set of the outputs of the flow-and-transport simulations. The 2D application cases refer to the highly non-stationary context of a contamination plume from a radioactive waste deposit located near the Chernobyl power station in Ukraine. A synthetic case is first examined and the estimation using real-world data is then discussed for comparing three estimation methods: usual kriging, external drift kriging with a simplified flow simulation as covariate, and kriging with numerical covariance.

Methods

"Numerical" covariance

In classical kriging approaches, the covariance is inferred from observations of the variable to be estimated. However kriging is more general and can be performed based on non-stationary covariances of type 𝐶(𝑥, 𝑥′). Especially, these covariances can be computed between each pair of grid nodes of the domain from a set of 𝑁 realizations of the variable [START_REF] Pannecoucke | Combining geostatistics and simulations of flow and transport to characterize contamination within the unsaturated zone[END_REF], such as:

𝐶(𝑥, 𝑥′) = 1 𝑁 ∑ [𝑍 𝑛 (𝑥) -𝑍(𝑥) ̅̅̅̅̅̅ ][𝑍 𝑛 (𝑥′) -𝑍(𝑥′) ̅̅̅̅̅̅̅ ] 𝑁 𝑛=1 (1)
where 𝐶(𝑥, 𝑥′) is the "numerical covariance" between grid nodes 𝑥 and 𝑥′ and 𝑍 𝑛 (𝑥) (respectively 𝑍 𝑛 (𝑥′)) is the value of the variable 𝑍 at grid node 𝑥 (respectively 𝑥′) for the 𝑛-th realization. For each grid node 𝑥, the "numerical" covariance at all 2D grid nodes 𝑥 ′ draws a 2D map. Fig. 1 shows the maps of numerical covariances related to two different points 𝑥, for a set of simulated concentration plumes (see further section 3.2.). It can be shown that the numerical covariance 𝐶(𝑥, 𝑥′) (or the numerical variogram 𝛾(𝑥, 𝑥′)) has the required mathematical properties of a (discrete) covariance or variogram function [START_REF] De Fouquet | Exercices corrigés de géostatistiqu[END_REF][START_REF] Pannecoucke | Combinaison de la géostatistique et des simulations à base physiqueapplication à la caractérisation de panaches de contaminants[END_REF].

For the application (Section 3) the realizations of the variable of interest are generated by flow-andtransport simulations of a contaminant plume. First, a set of geostatistical simulations of the hydraulic conductivity is constructed and simulations of the source term of the contamination are taken from a previous work [START_REF] Nguyen | Contribution à la quantification des incertitudes portées par la variabilité spatiale des déchets radioactifs enterrés à Tchernobyl[END_REF]. The flow-and-transport simulations are based on the synthesis of several studies (Le [START_REF] Coz | Do transient hydrological processes explain the variability of strontium-90 activity in groundwater downstream of a radioactive trench near Chernobyl?[END_REF]. Nevertheless, it should be noted that the influence of the range of the hydrogeological parameters was reduced compared with the variability between realizations and with the boundary conditions [START_REF] Pannecoucke | Impact of spatial variability in hydraulic parameters on plume migration within unsaturated surficial formations[END_REF][START_REF] Pannecoucke | Combinaison de la géostatistique et des simulations à base physiqueapplication à la caractérisation de panaches de contaminants[END_REF]).

The kriging with numerical covariance can be seen as a generalization of kriging in natural flow coordinates. Indeed, a stationary covariance 𝐾 applied after a space transformation 𝑓 can be written as 𝐾(𝑓(𝑥) -𝑓(𝑥 ′ )) and therefore is of the form 𝐶(𝑥, 𝑥') when the function 𝑓 is not linear. Conversely, in the case of a complex phenomenon, such as the development and movement of a plume, it is not necessarily the case that a space transformation 𝑓 and a stationary covariance 𝐾 exist such that the non-stationary covariance 𝐶(𝑥, 𝑥') is of the form 𝐾(𝑓(𝑥) -𝑓(𝑥 ′ )). It is thus pragmatic to calculate the empirical non-stationary covariance (the "numerical covariance") on a set of outputs of the numerical flow-and-transport simulations.

Non-stationary Gaussian transformation

In practice, kriging with a numerical covariance does not ensure that the estimated concentration is positive. The non-linear estimation in the well-known anamorphosed Gaussian model [START_REF] Chilès | Geostatistics: Modeling Spatial Uncertainty[END_REF] allows avoiding negative estimates for positive quantities.

Due to the non-stationarity of the local distributions of the simulated contaminant concentrations (especially, inside vs. outside the contaminant plume), a Gaussian transformation of these distributions is performed at each grid node. The non-stationary anamorphosis 𝛷 𝑥 is calculated from the 𝑁 outputs of the flow-and-transport simulation at the grid node 𝑥. Thus:

𝑍 𝑛 (𝑥) = 𝛷 𝑥 [𝑌 𝑛 (𝑥)] (2) and conversely 𝑌 𝑛 (𝑥) = 𝜙 𝑥 -1 [𝑍 𝑛 (𝑥)]
where 𝑌 𝑛 (𝑥) is the Gaussian transform of the contaminant concentration 𝑍 𝑛 (𝑥) of the 𝑛-th flow-andtransport simulation at grid node 𝑥. Here the anamorphosis is calculated using its Hermite polynomial expansion (with 30 polynomials; [START_REF] Chilès | Geostatistics: Modeling Spatial Uncertainty[END_REF].

The numerical covariance is then calculated on the set of the 𝑁 fields of Gaussian transformed variables.

The Gaussian data used for the estimation are deduced from the measurements (or from the reference data in the synthetic case) using the inverse of the anamorphosis (𝜙 𝑥 -1 ). The kriging estimate and the kriging variance are calculated in the Gaussian world.

Finally the conditional expectation of the concentration, its quantiles or any related functions, are calculated through the known relationship with the spatial distribution of a Gaussian random function; for example [START_REF] Matheron | A Simple Substitute for Conditional Expectation: The Disjunctive Kriging[END_REF][START_REF] Chilès | Geostatistics: Modeling Spatial Uncertainty[END_REF]:

𝐸(𝑍(𝑥)|𝑍(𝑥 𝑖 ), 𝑖 = 1, … , 𝑛) = ∫ 𝜙 𝑥 (𝑌 * (𝑥) + 𝜎(𝑥)𝑢) 𝑔(𝑢)𝑑𝑢 (3)
where 𝐸(𝑍(𝑥)|𝑍(𝑥 𝑖 ), 𝑖 = 1, … , 𝑛) is the conditional expectation of the concentration; 𝑌 * (𝑥) is the kriging with known mean (equal to 0), 𝜎(𝑥) the standard deviation of the kriging error, 𝑔 the density of the reduced normal distribution and u is the variable to be integrated. This expression can be calculated using the formula of the Hermite-polynomial expansion of the Gaussian anamorphosis [START_REF] Chilès | Geostatistics: Modeling Spatial Uncertainty[END_REF]. A known correction is applied on the standard deviation if the kriging with unknown mean (ordinary kriging) is used [START_REF] Emery | Ordinary multigaussian kriging for mapping conditional probabilities of soil properties[END_REF]. Equation (3) ensures that the estimated concentration is positive, because the concentrations outputs of the flow-and-transport simulations used for calculating the anamorphosis are positive.

Estimation workflow with numerical covariances

In summary, the method involves the following steps:

 simulating a set of unconditional hydraulic conductivity fields (geostatistical unconditional simulations by any classical method, e.g. the turning bands as used here; Chilès and Delfiner, 2012);

 simulating initial source term of the contamination (geostatistical conditional simulations in the following example)

 simulating flow-and-transport on each hydraulic conductivity field to obtain a set of concentration fields;

 calculating the non-stationary anamorphosis of the concentrations at each grid node 𝛷 𝑥 from the set of simulated concentration fields;

 calculating the Gaussian variable for each simulated concentration field using 𝛷 𝑥 -1 ;

 calculating the numerical covariances from the set of Gaussian variables;

 calculating Gaussian data as the Gaussian transform of concentrations data using 𝛷 𝑥 -1 at data points;

 kriging of Gaussian variable from Gaussian data and computing associated kriging variance, using the numerical covariances on Gaussian transforms;

 returning to the estimated concentration (or any function of the concentration, e.g. conditional quantiles or conditional distribution) thanks to Eq. ( 3)

Case study

The trench T22 in Chernobyl

The Chernobyl accident of April 1986 released large amounts of radionuclides into the surrounding environment. In 1986 and 1987, emergency clean-up was carried out in the irradiated Red Forest, in the vicinity of the Chernobyl nuclear power plant. About 500×10 3 m 3 of contaminated topsoil, vegetation and other radioactive materials were buried in situ in hundreds of trenches dug within the permeable sandy superficial formation, directly above the water table and without engineered moisture-proof barriers [START_REF] Kashparov | Radionuclide migration in the experimental polygon of the Red Forest waste site in the Chernobyl zone -Part 1: Characterization of the waste trench, fuel particle transformation processes in soils, biogenic fluxes and effects on biota[END_REF].

In 1998, an experimental research facility was built around the trench T22 in order to investigate radionuclide migration processes. Then, for about 20 years, the trench T22 has been monitored extensively, allowing an advanced characterization of the site geology, the water flow regime in the unsaturated and saturated zones, the radionuclide distribution, and the speciation within and outside the trench (e.g., [START_REF] Bugai | Radionuclide migration at experimental polygon at Red Forest waste site in Chernobyl zone. Part 2: Hydrogeological characterization and groundwater transport modeling[END_REF][START_REF] Kashparov | Radionuclide migration in the experimental polygon of the Red Forest waste site in the Chernobyl zone -Part 1: Characterization of the waste trench, fuel particle transformation processes in soils, biogenic fluxes and effects on biota[END_REF]. Moreover, acquired data and knowledge were integrated into numerical models that simulate the migration of radionuclides, especially 90 Sr (Strontium-90), from the trench to the underlying aquifer [START_REF] Bugai | Radionuclide migration at experimental polygon at Red Forest waste site in Chernobyl zone. Part 2: Hydrogeological characterization and groundwater transport modeling[END_REF][START_REF] Nguyen | Contribution à la quantification des incertitudes portées par la variabilité spatiale des déchets radioactifs enterrés à Tchernobyl[END_REF]. Le [START_REF] Coz | Do transient hydrological processes explain the variability of strontium-90 activity in groundwater downstream of a radioactive trench near Chernobyl?[END_REF] present a detailed synthesis of previous works used for the flow-and-transport modeling near the trench.

Flow-and-transport simulations

Between 2000 and 2015, groundwater was collected for 90 Sr activity measurements generally once or twice a year at four sampling depths (layers L1 to L4, from top to bottom) in six wells (AP1 to AP6, from upstream to downstream) along a cross-section A-B, approximatively perpendicular to the trench T22 and assumed to be parallel to the groundwater flow lines (Fig. 2 a). A groundwater flow-and-transport model was used for simulating in 2D, i.e. along the vertical cross-section A-B (Fig. 2 b), the development of the 90 Sr plume highlighted by this monitoring (Le [START_REF] Coz | Do transient hydrological processes explain the variability of strontium-90 activity in groundwater downstream of a radioactive trench near Chernobyl?[END_REF]. Model parameters were constant for each geological unit: eolian deposits, alluvial deposits with floodplain or channel facies (Le [START_REF] Coz | Do transient hydrological processes explain the variability of strontium-90 activity in groundwater downstream of a radioactive trench near Chernobyl?[END_REF]. These parameters, as well as the boundary conditions, were set such that the simulated time series of 90 Sr activity were consistent with most of the measurements. This simulation with especially uniform hydraulic conductivity by geological unit is called the calibrated simulation hereafter.

For this study, 400 new flow-and-transport simulations are run by considering the uncertainties on the model parameters (e.g., the parameter that governs radionuclide retention on aquifer material) and boundary conditions (e.g., the percolation rate through the trench T22). In addition, for each run,  the spatial variability of hydraulic conductivity is simulated through a random field for which the properties are determined based on textural analysis of aquifer material samples (eolian and alluvial sand). T

 the initial distribution of 90 Sr activity within the trench T22 is a geostatistical conditional simulation (Fig. 2 b) constructed from both in situ emission profiles and laboratory analysis of waste samples [START_REF] Nguyen | Contribution à la quantification des incertitudes portées par la variabilité spatiale des déchets radioactifs enterrés à Tchernobyl[END_REF].

Spatial estimation

A set of 357 simulated fields (some simulations did not converge) of the 90 Sr activity in the saturated zone along the cross-section A-B is extracted from the flow-and-transport simulations on 1 st October 2006 (01/10/2006). On this date, the 90 Sr activity was measured in groundwater collected at all sampling depths in wells AP1 to AP6 (Fig. 2 b). On the calibrated simulation, the 90 Sr plume extends up to 40 m from the trench with a downward component (Fig 3a). The geometry of the plume in the simulation set is similar, although its extent varies significantly (e.g., Fig 3b). The distribution of the 90 Sr activities computed from the simulation set at the sampling wells shows significant ranges, especially in the core of the plume (highest values) (Fig. 4), yet these ranges generally include the activity of the calibrated simulation and of the measurements.

Three estimation methods of the 90 Sr activity are compared in a cross-validation approach: (i) ordinary kriging (OK); (ii) kriging with an external drift (KED), the drift being the 90 Sr activity field extracted from the calibrated simulation (Fig. 3a); and (iii) kriging with numerical covariances (KNC) using a nonstationary Gaussian transformation, with covariances computed on the Gaussian transforms of the simulated fields of 90 Sr activity (simulations 2 to 357).

First, the simulated field 1 (Fig. 3b) is used as a synthetic case in order to compare the results with a known "reality" (this field is not considered for the covariance computation). Estimations are conditioned (i) at the wells AP1 to AP6 and then (ii) only at the most downstream well (AP6), the 90 Sr activity at the wells AP1 to AP5 being used for the assessment.

Then the actual measurements are considered. As for the synthetic case (ii) the estimations are performed only with the data belonging to the most downstream well (AP6), the 90 Sr activity at the wells AP1 to AP5 being used for the assessment.

Results

Synthetic case

The 90 Sr activities extracted from the simulated field 1 (i.e. the synthetic data) at the sampling locations are relatively high with respect to the activities extracted from the other simulations, included the calibrated one (Fig. 4). At the wells, the relationship between activity "data" and covariate (calibrated simulation) is approximately linear (Fig. 4).

Pragmatically, OK and KED kriging are done using an isotropic linear variogram, for avoiding the variogram fitting phase under unfavorable conditions (few data and high variability). When all data at the wells AP1 to AP6 are used in the estimation, the OK estimation (Fig 5a) is consistent with the simulated field 1 (Fig. 3b) only close to these wells, whereas the KED and the KNC estimations (Fig. 5 b andc) reproduce the geometry (extent and downward component) of the 90 Sr plume on the whole modeling domain. In addition, OK and KED result in negative 90 Sr activity estimated values in the bottom part of the domain, below the 90 Sr plume. KED and KNC are similar according to the root mean square error (RMSE) criterion (Tab. 1).

When only synthetic data at the well AP6 are used, the geometry of the 90 Sr plume is still well reproduced by KED and KNC but not at all by the OK (Fig. 6). KED results in negative concentrations outside the 90 Sr plume. The RMSE criterion calculated on all nodes (Tab. 1), or on the five validation wells only, favors the KNC, and it is nearly the same for the other statistical criteria (Fig. 7 and Tab 2).

A scatter diagram between the estimate and "reality" (Fig. 8) confirms the good accuracy of KED and KNC.

A test with the average of the 356 others simulations (instead of the calibrated simulation) as external drift leads to very similar results, only slightly less precise.

Actual measurements

The actual measured 90 Sr activities at some sampling locations differ significantly from the activities extracted from the calibrated simulation, but are generally (except in one location) included in the range of activities extracted from the 357 additional simulations (Fig. 4). The relationship between the activity at the wells and the covariate is no longer linear.

As for the synthetic case, when the estimations are conditioned by the actual measurements at the well AP6, OK does not reproduce the lateral variability of the 90 Sr activity with depth (Fig. 9a). KED maintains the geometry of the 90 Sr plume of the calibrated simulation, although the simulated activity is about a half order of magnitude lower (Fig. 9b). The geometry of the 90 Sr plume estimated by KNC differs significantly: its extent is reduced but the 90 Sr activity in its core is higher (Fig. 9c). At the five validation wells, all estimates are positive. The relative performances of KED and KNC are here more difficult to decide, because the "best" method varies according to the criterion (Fig. 10 and Tab. 3): the mean error is smaller for KNC, the median error and the RMSE are similar for the two methods, whereas the mean absolute error favors the KED. The scattered points between the estimates and observations are in this case further away from the first bisector (Fig. 11). This confirms the greater inaccuracy than for the synthetic case.

Discussion and conclusion

The main processes that govern the development of a contaminant plume in groundwater can be simulated using a flow-and-transport model. However, due to uncertainties in input data or misrepresentation of some sub-processes, model output cannot be considered as an accurate representation of the spatial distribution of the contaminant in groundwater at a given date. Usual kriging, e.g., ordinary kriging (OK), offers an alternative, but a high density of observations is necessary for providing realistic estimations. When observations are scarce, geostatistical approaches that take account of physically-based information, e.g., kriging with external drift (KED) or kriging with numerical covariance (KNC), are more relevant.

In KED, physically-based information consists of one or several drifts, which governs the geometry of the estimated contaminant plume. Following the literature, the output of a simplified flow model, namely with constant conductivity, is chosen here as external drift. KED requires, therefore, only one flow-and-transport simulation. In the case of actual observations, this drift is not necessarily consistent with those data; then KED could modify the variability of the estimated concentrations (location and level) of the plume.

With KNC, a set of several hundred flow-and-transport simulations are used as physically-based information for computing the numerical covariances (and optionally the non-stationary anamorphosis). On synthetic cases, when the observations are consistent with the simulation set (i.e., when observations are extracted from one reference simulation), the geometry of the contaminant plume of the reference simulation is fairly well reproduced. In the case of actual observations, the geometry of the contaminant plume could significantly differ from that of the simulated plumes. If the numerical covariance does not reflect this discrepancy, then the KNC will lead to inaccurate results.

For the actual measurements of the trench T22, the estimated concentrations seem more contrasted with KNC than with KED (Fig. 9), with larger maximum values. In addition, in all cases, concentration fields estimated by KNC associated with a Gaussian anamorphosis are strictly positive.

These results corroborate the conclusions of the previous work of [START_REF] Pannecoucke | Combining geostatistics and simulations of flow and transport to characterize contamination within the unsaturated zone[END_REF], which focused on a synthetic case that represents the development of a contaminant plume in the unsaturated zone. Indeed, this study, based on a real dataset from the saturated zone, confirms that the non-stationary numerical covariances allow integrating phenomenological knowledge in the kriging estimation. Moreover, it also shows that kriging with numerical covariance could be more accurate than external drift kriging, which does not consider the uncertainties related to this covariate, especially when the geometry of the contaminant plume is complex. The quality of the estimation depends on the ability of the physically-based realizations used for computing the numerical covariances to be consistent with the actual flow-and-transport processes and the corresponding uncertainties.

However, the (unique) simplified flow-and-transport simulation used as drift for the KED requires a much lower calculation time than the set of simulations needed for computing the numerical covariances. One reviewer points out that "the workflow needed for the KNC also requires a consistent number of additional steps as compared with the other alternatives."

Rather than a new "better method", kriging with numerical covariance appears to be an alternative to external drift kriging, which enriches the estimation methods useful to the practitioner. This general method can be applied to diverse issues in hydrogeology, in multivariate or spatio-temporal modeling, when estimating or conditioning simulations of hydraulic head or fluxes, backtracking of a plume to locate the source of a pollution [START_REF] Pannecoucke | Combinaison de la géostatistique et des simulations à base physiqueapplication à la caractérisation de panaches de contaminants[END_REF], and in hydro-geochemical modelling. The method obviously applies to other contexts (e.g. air quality), and provides a very general class of covariances, e.g. on graphs or other non-Euclidean geometry. Figure 7. Synthetic case. Profiles of estimated activity by OK, EDK or KNC at wells AP1 to AP5, conditioned by "data" at well AP6 and "reality". Logarithmic scale.

Figure_8. Synthetic case. Scatter diagram between estimate (abscissa) and reference value (ordinate) for the three estimation methods (OK, KED, KNC) at the five validation wells for the estimation from AP6.

  Figure captionsFigure 1a,b. Examples of numerical covariance maps, for two different points 𝑥 (red circles).

Figure 2 .

 2 Figure 2. Views of the 3-D model of the 90 Sr inventory in 1987 in the trench T22 (computed from the mean inventory estimated by Nguyen, 2017) and locations of the wells monitored for groundwater activity. (a) Horizontal view with cumulative 90 Sr activity in depth; (b) vertical cross-section A-B .

Figure 3 .

 3 Figure 3. Activity of 90 Sr in the saturated zone downstream of the cross-section A-B of the trench T22 simulated with the flow-and-transport model on 01/10/2006: (a) calibrated simulation used as external drift; and (b) simulation 1 taken as reference for the synthetic case study Figure 4. 90 Sr activity at the sampling locations. Scatter diagram between (a) the synthetic case in abscissa (horizontal axis) and the calibrated simulation in ordinate (vertical axis) and (b) the actual measurements (horizontal axis) and the calibrated simulation in ordinate, . The grey bar represents the range of the set of 357 simulations. The black line corresponds to the first bisector.

Figure 5 .

 5 Figure5. Synthetic case (i.e., flow-and-transport simulation 1). Activity of 90 Sr in the saturated zone downstream of the cross-section A-B of the trench T22 (Fig.2) estimated on 01/10/2006 conditioned at wells AP1 to AP6: (a) ordinary kriging (OK); (b) kriging with external drift (KED); and (c) kriging with numerical covariances using a non-stationary Gaussian transformation (KNC). The non-filled areas (in a and b) correspond to negative estimated activities.

Figure 6 .

 6 Figure 6. Synthetic case. Activity of 90 Sr in the saturated zone downstream of the cross-section A-B of the trench T22 (Fig. 2) estimated on 01/10/2006 conditioned at well AP6 only: (a) ordinary kriging (OK); (b) kriging with external drift (KED); and (c) kriging with numerical covariances using a nonstationary Gaussian transformation (KNC). The non-filled area (in b) corresponds to negative estimated activities.

Figure 9 .

 9 Figure 9. Activity of 90 Sr in the saturated zone downstream of the cross-section A-B of the trench T22 (Fig. 2) estimated on 01/10/2006 conditioned with measurements at well AP6: (a) ordinary kriging (OK); (b) kriging with external drift (KED); and (c) kriging with numerical covariances using a nonstationary Gaussian transformation (KNC). The non-filled area (in b) corresponds to negative estimated activities.

Figure 10 .

 10 Figure 10. Profiles of estimated activity by OK, EDK or KNC at wells AP1 to AP5, conditioned by actual measurements at well AP6, and validation measurements. Logarithmic scale. Figure 11. Actual measurements. Scatter diagram between estimate (abscissa) and exact value (ordinate) for the three estimation methods (OK, KED, KNC) at the five validation wells for the estimation from AP6.

Figure 1

 1 Figure 1

Figure 2

 2 Figure 2 Click here to access/download;Figure;Fig2R.pdf

Figure 3

 3 Figure 3 Click here to access/download;Figure;Fig3R.pdf

Figure 4

 4 Figure 4Click here to access/download;Figure;Fig4R.pdf

Figure 5

 5 Figure 5 Click here to access/download;Figure;Fig5R.pdf

Figure 6

 6 Figure 6Click here to access/download;Figure;Fig6R.pdf

Figure 7

 7 Figure 7

Figure 8

 8 Figure 8

Figure 9

 9 Figure 9Click here to access/download;Figure;Fig9R.pdf

Figure 10

 10 Figure 10

Figure 11

 11 Figure 11

Table 1 .

 1 RMSE at all 1785 grid nodes, for the two conditioning sets of wells. The very weakly negative values of KED are retained. The numbers have been rounded

	Tables
	Well set	OK	KED KNC
	6 conditioning wells	363. 102. 106.
	1 conditioning well (AP6) 658. 135. 94.

Table 2 .

 2 Cross-validation statistics at the five validation wells (AP1 to AP5, 20 values) estimated from the data of the AP6 well, for the synthetic case. Negative estimates are zeroed for the statistics of the error, but their number is reported. The error is calculated as estimate -exact data.

	Statistic	OK	KED(0)	KNC
	Mean error	-127.	-103.	3.
	Median error	306.	0.	10.
	Mean absolute error	692.	134.	87.
	RMSE	908.	227.	135.
	Number of negative estimates	0	9	0

Table 3 .

 3 Cross-validation statistics at the five validation wells (AP1 to AP5, 20 values) estimated from the data of the AP6 well, for the actual measurements.

	Statistic	OK	KED	KNC
	Mean error	-199.	-195.	-80.
	Median error	37.	0.	-3.
	Mean absolute error	312.	227.	267.
	RMSE	589.	472.	477.
	Number of negative estimates	0	0	0

Acknowledgments

The authors cordially thank the reviewers and editor whose comments greatly helped to improve the paper. Following the remarks of one reviewer, the "calibrated" simulation was used as external drift, in preference to the average of the set of simulations.

Funding Information and Conflicts of Interest

This study benefited from the Kri-Terres project, supported by the French National Radioactive Waste Management Agency (ANDRA) under the French "programme d'investissements d'avenir".

The authors declare there is no conflict of interest in connection with the submitted material


Valentina Baldazzi
email: valentina.baldazzi@inria.fr

Delphine Ropers

Jean-Luc Gouzé

Tomas Gedeon

Hidde De Jong
email: hidde.de-jong@inria.fr

Resource allocation accounts for the large variability of rate-yield phenotypes across bacterial strains

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

Introduction

Microbial growth consists of the conversion of nutrients from the environment into biomass. This flux of material is coupled with a flux of energy from the substrate to small energy cofactors (ATP, NADH, NADPH, ...) driving biomass synthesis forward and releasing energy in the process [START_REF] Schachter | Microbe[END_REF]. The growth of microorganisms has been profitably analyzed from the perspective of resource allocation, that is, the assignment of limiting cellular resources to the different 1 biochemical processes underlying growth [START_REF] Scott | Interdependence of cell growth and gene expression: Origins and consequences[END_REF][START_REF] Scott | Emergence of robust growth laws from optimal regulation of ribosome synthesis[END_REF]4,[START_REF] Giordano | Dynamical allocation of cellular resources as an optimal control problem: Novel insights into microbial growth strategies[END_REF][START_REF] Weiße | Mechanistic links between cellular tradeoffs, gene expression, and growth[END_REF][START_REF] Reimers | Cellular trade-offs and optimal resource allocation during cyanobacterial diurnal growth[END_REF][START_REF] Bosdriesz | How fast-growing bacteria robustly tune their ribosome concentration to approximate growth-rate maximization[END_REF][START_REF] Towbin | Optimality and suboptimality in a bacterial growth law[END_REF][START_REF] Maitra | Bacterial growth laws reflect the evolutionary importance of energy efficiency[END_REF][START_REF] Dourado | An analytical theory of balanced cellular growth[END_REF][START_REF] Metzl-Raz | Principles of cellular resource allocation revealed by condition-dependent proteome profiling[END_REF]. It is often considered that proteins, the main component of biomass, are also the bottleneck resource for growth. Proteins function as enzymes in carbon and energy metabolism and they constitute the molecular machines responsible for the synthesis of macromolecules, in particular proteins themselves. The composition of the proteome in a given growth condition can therefore be interpreted as the resource allocation strategy adopted by the cells to exploit available nutrients.

Two macroscopic criteria for characterizing microbial growth are growth rate and growth yield. The former refers to the rate of conversion of substrate into biomass, and the latter to the efficiency of the process, that is, the fraction of substrate converted into biomass. Several empirical relations between proteome composition vs growth rate and growth yield have been established. A linear relation between growth rate and the ribosomal protein fraction of the proteome holds over a large range of growth rates and for a variety of microbial species [START_REF] Scott | Interdependence of cell growth and gene expression: Origins and consequences[END_REF][START_REF] Forchhammer | Growth rate of polypeptide chains as a function of the cell growth rate in a mutant of Escherichia coli 15[END_REF][START_REF] Bremer | Modulation of chemical composition and other parameters of the cell by growth rate[END_REF]. Variants of this so-called growth law have been found for cases of reduced translation capacities [START_REF] Scott | Interdependence of cell growth and gene expression: Origins and consequences[END_REF] or different temperatures [START_REF] Herendeen | Levels of major proteins of Escherichia coli during growth at different temperatures[END_REF][START_REF] Mairet | Optimal proteome allocation and the temperature dependence of microbial growth laws[END_REF]. While the ribosomal protein fraction increases with the growth rate, the proteome fraction allocated to energy metabolism decreases [START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF][START_REF] Schmidt | The quantitative and condition-dependent Escherichia coli proteome[END_REF]. Moreover, within this decreasing fraction, Escherichia coli and other microorganisms move resources from respiration to fermentation pathways [START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF]. Simple mathematical models have been proposed to account for the above relations in terms of the requirements of self-replication of the proteome and the relative protein costs and ATP yields of respiration and fermentation [START_REF] Scott | Interdependence of cell growth and gene expression: Origins and consequences[END_REF]4,[START_REF] Giordano | Dynamical allocation of cellular resources as an optimal control problem: Novel insights into microbial growth strategies[END_REF][START_REF] Weiße | Mechanistic links between cellular tradeoffs, gene expression, and growth[END_REF][START_REF] Bosdriesz | How fast-growing bacteria robustly tune their ribosome concentration to approximate growth-rate maximization[END_REF][START_REF] Dourado | An analytical theory of balanced cellular growth[END_REF][START_REF] Mairet | Optimal proteome allocation and the temperature dependence of microbial growth laws[END_REF][START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF][START_REF] Mori | A yield-cost tradeoff governs Escherichia coli's decision between fermentation and respiration in carbon-limited growth[END_REF].

Most of these relations have been studied in experiments in which the same strain exhibits a range of growth rates in different environments, with different carbon sources. Even for a fixed environment, however, different strains of the same species may grow at very different rates and yields. For example, in a comparative study of seven E. coli strains, growth rates ranging from 0.61 to 0.97 h -1 , and (carbon) growth yields between 0.52 and 0.66, were observed during aerobic growth on glucose [START_REF] Monk | Multi-omics quantification of species variation of Escherichia coli links molecular features with strain phenotypes[END_REF]. Since the genes encoding enzymes in central carbon and energy metabolism are largely shared across the strains [START_REF] Monk | Multi-omics quantification of species variation of Escherichia coli links molecular features with strain phenotypes[END_REF], the yield differences are not due to different metabolic capacities but rather to different regulatory strategies, that is, different usages of the metabolic pathways of the cell. As another example, evolution experiments with E. coli have given rise to evolved strains that grow more than 40% faster, sometimes with higher growth yields, than the ancestor strain in the same environment [START_REF] Lacroix | Use of adaptive laboratory evolution to discover key mutations enabling rapid growth of Escherichia coli K-12 MG1655 on glucose minimal medium[END_REF]. Analysis of the underlying mutations reveals that the higher rates and yields of the evolved strains are not due to new metabolic capacities, but rather to modified regulatory strategies [START_REF] Lacroix | Use of adaptive laboratory evolution to discover key mutations enabling rapid growth of Escherichia coli K-12 MG1655 on glucose minimal medium[END_REF][START_REF] Utrilla | Global rebalancing of cellular resources by pleiotropic point mutations illustrates a multi-scale mechanism of adaptive evolution[END_REF].

Can the large variability of rate-yield phenotypes observed across different strains of the same species be explained by different resource allocation strategies, that is, different compositions of the proteome? In order to answer this question, we developed a coarse-grained resource allocation model that couples the fluxes of carbon and energy underlying microbial growth. The model was calibrated by means of existing data in the literature, without any parameter fitting, and its predictions were compared with a database of several hundreds of pairs of rates and yields of E. coli strains reported in the literature. The database includes wild-type strains as well as mutant strains obtained through directed mutagenesis or adaptive laboratory evolution (ALE).

We found that the predicted rates and yields in different growth conditions correspond very well with the observed variability of rate-yield phenotypes. The model also accounts for 2 the associated substrate uptake and acetate secretion rates. Contrary to what was expected from studies with a single strain in different conditions [START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF][START_REF] Zavřel | Quantitative insights into the cyanobacterial cell economy[END_REF][START_REF] Nanchen | Nonlinear dependency of intracellular fluxes on growth rate in miniaturized continuous cultures of Escherichia coli[END_REF][START_REF] Peebo | Proteome reallocation in Escherichia coli with increasing specific growth rate[END_REF], the model predicted no correlation between growth rate and acetate secretion rate in the case of different strains growing in the same environment, a result confirmed by the experimental data. A mapping from resource allocation to rate-yield phenotypes consistent with available proteomics data, however, could only be found when taking into account enzyme activities in addition to enzyme concentrations.

These results are interesting for both fundamental reseach and biotechnological applications. They show that the application of coarse-grained models can be used to predict multivariate phenotypes, without making any assumptions on optimality criteria, and reveal unexpected relations confirmed by the experimental data. The model is capable of predicting quantitative bounds on growth rates and yields within a specific environment, which can be exploited for rapidly screening performance limits of strains developed in synthetic biology and metabolic engineering.

Results

Coarse-grained model with coupled carbon and energy fluxes

Coarse-grained resource allocation models describe microbial growth by means of a limited number of macroreactions converting nutrients from the environment into proteins and other macromolecules. Several such models have been proposed, usually focusing on carbon fluxes [START_REF] Scott | Interdependence of cell growth and gene expression: Origins and consequences[END_REF]4,[START_REF] Giordano | Dynamical allocation of cellular resources as an optimal control problem: Novel insights into microbial growth strategies[END_REF][START_REF] Weiße | Mechanistic links between cellular tradeoffs, gene expression, and growth[END_REF][START_REF] Bosdriesz | How fast-growing bacteria robustly tune their ribosome concentration to approximate growth-rate maximization[END_REF][START_REF] Towbin | Optimality and suboptimality in a bacterial growth law[END_REF][START_REF] Mairet | Optimal proteome allocation and the temperature dependence of microbial growth laws[END_REF]. The coupling of carbon and energy fluxes is essential for understanding the relation between growth rate and growth yield. Few models, however, have taken into account the energy balance of microbial growth, that is, the use of substrate for the production of ATP and other energy cofactors fueling macromolecular synthesis and cellular maintenance. Among the notable exceptions, we cite the model of Basan et al. [START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF][START_REF] Mori | A yield-cost tradeoff governs Escherichia coli's decision between fermentation and respiration in carbon-limited growth[END_REF], which couples carbon and energy fluxes while abstracting from the reaction kinetics, and the model of Zavřel et al. [START_REF] Zavřel | Quantitative insights into the cyanobacterial cell economy[END_REF], which does provide such a kinetic view while restricting biomass to proteins.

Figure 1 presents a coarse-grained kinetic model that takes inspiration from and generalizes upon this previous work. While the model is generic, it has been instantiated for aerobic growth of E. coli in minimal medium with glucose or glycerol as the limiting carbon source. The model variables comprise cellular concentrations of proteins (p) and other macromolecules (DNA, RNA, lipids, and polysaccharides) (u) as well as central carbon metabolites (c) and ATP (a * ). All concentrations have units Cmmol gDW -1 , except for ATP [mmol gDW -1 ]. Five macroreactions are responsible for carbohydrate uptake and metabolism, ATP production by aerobic respiration and fermentation, and the synthesis of proteins and other macromolecules. The rates of the reactions, denoted by v mc , v mer , v mef , v r , and v q [Cmmol gDW -1 h -1 ], respectively, are defined by kinetic expressions involving protein, precursor metabolite, and energy cofactor concentrations. Details of the rate equations and the derivation of the model from basic assumptions on microbial growth can be found in Text S1.

The carbon entering the cell is included in the different biomass components or released in the form of CO 2 and acetate. CO 2 is produced by respiration and macromolecular synthesis, while acetate overflow is due to aerobic fermentation [START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF][START_REF] Gottschalk | Bacterial Metabolism[END_REF]. The carbon balance also includes

v r , v q , v mc , v mer , v mef , v d [Cmmol or mmol gDW -1 h -1 ],
the variables p, r, q, m c , m er , m ef , c, u, a * [Cmmol or mmol gDW -1 ], the resource allocation parameters φ q , φ r , φ c , φ er , φ ef [dimensionless], the degradation rate constant γ [h -1 ], the biomass density β [Cmmol gDW -1 ], the ATP yield and cost factors n mer , n mef , n r , n q [mmol Cmmol -1 ], and the CO 2 loss factors ρ mef , ρ rq [dimensionless] are formally defined in Text S1. The values of the parameters are derived in Text S2.

4 the turnover of macromolecules, which is responsible for a large part of cellular maintenance costs ( [START_REF] Van Bodegom | Microbial maintenance: a critical review on its quantification[END_REF] and Text S2).

The energy balance is expressed in terms of the production and consumption of ATP. While energy metabolism also involves other energy cofactors (NADP, NADPH, . . .), the latter can be converted into ATP during aerobic growth [START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF][START_REF] Gottschalk | Bacterial Metabolism[END_REF]. The ATP yields of respiration and fermentation (n mer and n mef ) as well as the ATP costs of the synthesis of proteins and other macromolecules (n r and n q ) are determined by the stoichiometry of the underlying metabolic pathways and the biomass composition ( [START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF][START_REF] Kaleta | Metabolic costs of amino and protein production in Escherichia coli[END_REF] and Text S2). When total ATP production and consumption in growing microbial cells are computed from n mer v mer + n mef v mef and n r v r + n q v q , respectively, the former usually largely exceeds the latter [START_REF] Feist | A genome-scale metabolic reconstruction for Escherichia coli K-12 MG1655 that accounts for 1260 ORFs and thermodynamic information[END_REF][START_REF] Russell | Energetics of bacterial growth: balance of anabolic and catabolic reactions[END_REF]. This so-called uncoupling phenomenon is explicitly accounted for by an energy dissipation term v d in the energy balance (Text S1).

Like in other resource allocation models, the proteome is subdivided into categories [START_REF] Scott | Interdependence of cell growth and gene expression: Origins and consequences[END_REF][START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF]. We distinguish ribosomes and other translation-affiliated proteins, enzymes in central carbon metabolism, enzymes in respiration and fermentation metabolism, and other (housekeeping) proteins, with concentrations r, m c , m er , m ef , and q, respectively. Each category of protein catalyzes a different macroreaction in Figure 1: ribosomes are responsible for protein synthesis, enzymes for carbon and energy metabolism, and housekeeping proteins for the synthesis of macromolecules other than proteins. The protein synthesis capacity of the cell, given by the total protein synthesis rate v r , is distributed over the protein categories using five fractional resource allocation parameters that sum to 1: φ q , φ r , φ c , φ er , and φ ef . Fixing the resource allocation parameters determines the model dynamics and therefore the growth phenotype [START_REF] Dourado | An analytical theory of balanced cellular growth[END_REF][START_REF] Zavřel | Quantitative insights into the cyanobacterial cell economy[END_REF][START_REF] De Groot | Elementary Growth Modes provide a molecular description of cellular self-fabrication[END_REF].

Contrary to most models of microbial growth, the biomass includes other cellular components (DNA, RNA, metabolites, . . .) in addition to proteins. The growth rate µ [h -1 ] directly follows from the biomass definition, under the assumption that the total biomass concentration 1/β is constant ( [START_REF] Jong | Mathematical modelling of microbes: metabolism, gene expression and growth[END_REF] and Text S1). The growth rate captures the specific accumulation of biomass corrected for maintenance:

µ = β (v mc -v mer -ρ mef v mef -(ρ rq -1) (v r + v q )) -γ, (1)
where ρ mef and ρ rq -1 denote the fractional loss of CO 2 by fermentation and macromolecular synthesis, respectively. The growth yield is defined as the fraction of carbon taken up by the cells and included in the biomass:

Y = 1 β µ v mc . ( 2
)
Yields are dimensionless and vary between 0 and 1. The model in Figure 1 was calibrated using data from the literature for batch or continuous growth of E. coli in minimal medium with glucose or glycerol. In brief, for the E. coli reference strain BW25113, we collected for each growth medium the growth rate and metabolite uptake and secretion rates [START_REF] Peebo | Proteome reallocation in Escherichia coli with increasing specific growth rate[END_REF][START_REF] Van Rijsewijk | Large-scale L13C-flux analysis reveals distinct transcriptional control of respiratory and fermentative metabolism in Escherichia coli[END_REF][START_REF] Gerosa | Pseudo-transition analysis identifies the key regulators of dynamic metabolic adaptations from steady-state data[END_REF], as well as protein and metabolite concentrations [START_REF] Schmidt | The quantitative and condition-dependent Escherichia coli proteome[END_REF][START_REF] Gerosa | Pseudo-transition analysis identifies the key regulators of dynamic metabolic adaptations from steady-state data[END_REF][START_REF] Bennett | Absolute metabolite concentrations and implied enzyme active site occupancy in Escherichia coli[END_REF]. Using additional assumptions based on literature data, we fixed a unique set of parameters for each condition (batch vs. continuous growth, glucose vs. glycerol), without parameter fitting (Text S2). The resulting set of quantitative models provides a concise but comprehensive representation of microbial growth in different environments.

Predicted and observed rate-yield phenotypes for Escherichia coli

Using the model in Figure 1, we explored the effect of alternative resource allocation strategies on growth rate and growth yield. A resource allocation strategy is defined by the values for the resource allocation parameters: (φ q , φ r , φ c , φ er , φ ef ) (Figure 2A). The medium-specific instantiations of the model were used to determine how the rate and yield change, within a given condition, when resource allocation is allowed to vary over the entire range of possible strategies. We therefore uniformly sampled the space of possible combinations of resource allocation parameters, while keeping the fraction attributed to housekeeping proteins φ q constant [START_REF] Scott | Interdependence of cell growth and gene expression: Origins and consequences[END_REF]. The model was considered at steady state, corresponding to balanced growth of the culture. The steady state for each resource allocation strategy was determined numerically (Methods).

Figure 2B-C shows the space of possible rate-yield phenotypes for batch growth on glucose or glycerol. A first observation is that the predicted combinations of rate and yield are bounded. The growth rate does not exceed 1.1 h -1 , and for all but the lowest growth rates, the growth yield is larger than 0.3. A second observation is that, for low growth rates, the maximum growth yield increases with the rate, whereas it decreases for high growth rates. The initial yield increase is due to the proportionally lower burden of the maintenance costs at higher growth rates [START_REF] Pirt | The maintenance energy of bacteria in growing cultures[END_REF]. The subsequent decrease of the maximum yield reflects a trade-off that has been much investigated in microbial physiology and ecology [START_REF] Lipson | The complex relationship between microbial growth rate and yield and its implications for ecosystem processes[END_REF][START_REF] Beardmore | Metabolic trade-offs and the maintenance of the fittest and the flattest[END_REF].

The model predictions can be directly compared with experimental data. We compiled a database of measured rates and yields reported in the literature for growth of E. coli in minimal medium with glucose or glycerol (Files S1 and S2), and plotted the measurements in the predicted phenotype spaces (Figure 2B-C). The database includes the reference wildtype strain used for model calibration (BW25113), other E. coli wild-type strains, strains with mutants in regulatory genes, and strains obtained from ALE experiments. Apart from the rate and yield of the reference strain [START_REF] Van Rijsewijk | Large-scale L13C-flux analysis reveals distinct transcriptional control of respiratory and fermentative metabolism in Escherichia coli[END_REF], none of the data points plotted in Figure 2B-C were used for calibration.

The measured rates and yields correspond very well with the model predictions: almost all data points fall inside the predicted phenotype space and much of the space is covered by the data points. The range of high growth rates is enriched in data points, which reflects the bias that E. coli wild-type and mutant strains grow relatively fast on glucose and glycerol, and that in most ALE experiments the selection pressure is tilted towards growth rate. Interestingly, the highest growth rates on glucose attained in ALE experiments, just above 1 h -1 [START_REF] Lacroix | Use of adaptive laboratory evolution to discover key mutations enabling rapid growth of Escherichia coli K-12 MG1655 on glucose minimal medium[END_REF][START_REF] Monk | iML1515, a knowledgebase that computes Escherichia coli traits[END_REF], agree quite well with the highest predicted growth rates.

The BW25113 strain used for calibration has a low growth yield on glucose (equal to 0.50, [START_REF] Van Rijsewijk | Large-scale L13C-flux analysis reveals distinct transcriptional control of respiratory and fermentative metabolism in Escherichia coli[END_REF]). Many deletion mutants of regulatory genes somewhat increase that yield [START_REF] Van Rijsewijk | Large-scale L13C-flux analysis reveals distinct transcriptional control of respiratory and fermentative metabolism in Escherichia coli[END_REF], but still fall well below the maximally predicted yield. The growth yield of some other wild-type strains is significantly higher, for example the W strain achieves a yield of 0.66 at a growth rate of 0.97 h -1 [START_REF] Monk | Multi-omics quantification of species variation of Escherichia coli links molecular features with strain phenotypes[END_REF]. The highest growth yield is achieved by an evolved strain (0.81, [START_REF] Zamboni | Multidimensional optimality of microbial metabolism[END_REF]), agreeing quite well with the maximum predicted growth yield for that growth rate. The latter strain does not secrete any acetate while growing on glucose [START_REF] Zamboni | Multidimensional optimality of microbial metabolism[END_REF], which contributes to the higher yield.

Similar observations can be made for growth of E. coli on glycerol, although in this case less experimental data points are available. The model predicts that the highest growth rate on glycerol is close to the highest growth rate on glucose, which is confirmed by experimental 6 ). The predicted yields are shown for the dilution rates ±10%. All predictions were made using the model in Figure 1, calibrated for the different growth conditions, and varying the resource allocation parameters between 0 and 1, except φ q (40,000-90,000 samples). The measurements of rate and yield reported in the source literature have been converted to units h -1 (growth rate) and a dimensionless unit corresponding to Cmmol biomass Cmmol -1 substrate (growth yield) (see Text S1 and Files S1 and S2 for details).

data [START_REF] Andersen | Are growth rates of escherichia coli in batch cultures limited by respiration?[END_REF]. In addition to batch growth, we also considered continuous growth in a chemostat. This required a recalibration of the model, since the environment is not the same as for batch growth (Text S2). Figure 2D shows the predicted rate-yield phenotype space for dilution rates around 0.2, 0.35, and 0.5 h -1 , as well as the observed rates and yields. Again, there is good correspondence between model and data. Most chemostat experiments reported in the literature have been carried out with the BW25113 and MG1655 wild-type strains. This absence of mutants and evolved strains may lead to an underestimation of the range of observed growth yields.

In the above comparisons of the model with the data, we made the assumption that the strains considered have the same metabolic capacities as the reference strain. That is, they utilize the same pathways in central carbon and energy metabolism to convert substrate into precursor metabolites and ATP. This assumption was satisfied by restricting the database to wild-type strains with essentially the same central carbon and energy metabolism [START_REF] Monk | Multi-omics quantification of species variation of Escherichia coli links molecular features with strain phenotypes[END_REF], mutant strains with deletions of genes encoding regulators instead of enzymes [START_REF] Van Rijsewijk | Large-scale L13C-flux analysis reveals distinct transcriptional control of respiratory and fermentative metabolism in Escherichia coli[END_REF], and short-term ALE mutants which have not had the time to develop new metabolic capacities [START_REF] Monk | iML1515, a knowledgebase that computes Escherichia coli traits[END_REF]. We also made the assumption that the parameter values are the same for all strains, so that differences in resource allocation strategies are the only explanatory variable. It is remarkably that, despite this strong assumption, the model of Figure 1 accounts very well for the observed variability of rate-yield phenotypes in E. coli.

Predicted and observed flux phenotypes for Escherichia coli

Growth rate and growth yield are defined in terms of carbon and energy fluxes through the population (Eqs. 1 and 2). Like rate and yield, some of these fluxes, in particular uptake and secretion rates, have been found to largely vary across E. coli strains growing in minimal medium with glucose [START_REF] Monk | Multi-omics quantification of species variation of Escherichia coli links molecular features with strain phenotypes[END_REF][START_REF] Lacroix | Use of adaptive laboratory evolution to discover key mutations enabling rapid growth of Escherichia coli K-12 MG1655 on glucose minimal medium[END_REF]. Can our model also reproduce the observed flux variability? We compared the model predictions of the uptake and secretion rates with their measurements, when available (File S1).

Figure 3A-B shows the projection of the predicted and measured glucose uptake rates on the rate-yield phenotype space. The correspondence of the model and the data is very high: more than 90% of the variance is captured by the model (Figure S1). This high correlation is not surprising. For a given growth rate µ and growth yield Y , the glucose uptake rate v mc is determined by the relation in Eq. 2. What the correspondence shows is that the total biomass concentration 1/β does not vary much across the different strains, agreeing with what was observed for steady-state growth of a given E. coli strain in different media [42,[START_REF] Basan | Inflating bacterial cells by increased protein synthesis[END_REF].

Another observable flux phenotype is the acetate secretion rate, which is an indicator of the functioning of energy metabolism. In aerobic conditions, E. coli has two different modes of ATP production: respiration and fermentation. Glucose and glycerol are taken up by the cells and degraded in the glycolysis pathway, eventually producing acetyl-CoA. Whereas acetyl-CoA enters the tricarboxylic acid (TCA) cycle in the case of respiration, it is secreted in the form of acetate during fermentation. In both cases, NADP and other reduced compounds are produced along the way and their recycling is coupled with the generation of a proton gradient across the membrane, enabling the production of ATP. Respiration is the more efficient of the two ATP production modes: in E. coli, respiration yields 26 ATP molecules per molecule of glucose and fermentation only 12 [START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF]. Glucose uptake rate (v mc ) for every predicted rate-yield combination in the case of batch growth of E. coli on minimal medium with glucose (Figure 2B). B. Idem for measured combinations of glucose uptake rate, growth rate, and growth yield (File S1). The predicted and observed secretion rates are strongly correlated (R 2 = 0.93, Figure S1), in agreement with Eq. 2. C. Relative acetate secretion rate (v mef /v mc ) for every predicted rate-yield combination. D. Idem for measured combinations of glucose uptake rate, acetate secretion rate, growth rate, and growth yield (File S1). The predicted and observed relative acetate secretion rates are correlated (R 2 = 0.66, Figure S1). There is no evidence in the data of a correlation between growth rate and relative acetate secretion rate (R 2 = 0.08, Figure S1). In addition to the relative acetate overflow, the model also captures the absolute acetate overflow very well (Figure S2).

9

How does the contribution of fermentation to ATP production change over the rate-yield phenotypes? Contrary to the glucose uptake rate, the growth rate and growth yield do not uniquely determine the acetate overflow rate v mef (Eqs. 1 and 2). Instead, the acetate overflow rate is the consequence of the resource allocation strategy adopted by the cell. We examined the predicted and observed contributions of fermentation to ATP production in the case of batch growth on glucose. In particular, we quantified the extent of fermentation by the relative rate of acetate overflow, that is, the fraction of carbon taken up that is secreted as acetate (v mef /v mc ). As shown in Figure 3C, over almost the entire range of growth rates, a continuum of ATP production, from pure respiration to combinations of respiration and fermentation, is possible according to the model. However, the ATP production mode is directly correlated with the growth yield. Maximum yield requires respiration without fermentation, whereas minimum yield is attained for maximum fermentation, where more than 30% of the carbon entering the cell is lost due to acetate overflow.

The predicted pattern of the variation of acetate overflow over the rate-yield phenotype space is largely confirmed by the experimental data (Figure 3D). Some of the fastest growing E. coli wild-type strains have no acetate overflow, like the W strain [START_REF] Monk | Multi-omics quantification of species variation of Escherichia coli links molecular features with strain phenotypes[END_REF], while some of the evolved strains grow very fast but with little acetate overflow as compared to their ancestors [START_REF] Zamboni | Multidimensional optimality of microbial metabolism[END_REF]. Respiration phenotypes are over-represented at high growth yields, over the whole range of growth rates, whereas fermentation phenotypes are associated with low growth yields. The model accounts for maximally 66% of the variance of the observed relative acetate overflow (Figure S1). This correspondence may be underestimated due to a bias in our definition of fermentation phenotypes, which is based on acetate overflow only and does not take into account the possible secretion of other by-products (formate, succinate, lactate, ...). Overflow of other fermentation products is known to occur during batch growth on glucose, especially for mutants [START_REF] Wolfe | The acetate switch[END_REF][START_REF] Enjalbert | Physiological and molecular timing of the glucose to acetate transition in Escherichia coli[END_REF][START_REF] Waegeman | Effect of iclR and arcA knockouts on biomass formation and metabolic fluxes in Escherichia coli K12 and its implications on understanding the metabolism of Escherichia coli BL21 (DE3)[END_REF].

Given the higher ATP yield of respiration, it is not surprising that the highest growth yields are obtained when respiration is preferred to fermentation. What might not have been expected, however, is that resource allocation strategies enabling cells to bypass fermentation at high growth rates are possible at all. It is well-known that, when growing an E. coli strain in different environments supporting increasingly higher growth rates, the contribution of fermentation to ATP production increases at the expense of respiration, as witnessed by the relative increase of acetate secretion ( [START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF] and Figure S3). In the scenario considered here, however, both the model and the data suggest that there exist strategies enabling strains to attain growth rates close to the maximum without shifting resources from respiration to fermentation.

Differences in enzyme activity are necessary to explain variations in protein abundance

The simple model of Figure 1 accounts remarkably well for the variability of observed rateyield and flux phenotypes in E. coli. It raises the question how the phenotypes are related to the underlying resource allocation strategies. This relation is not straightforward due to the multiple feedback loops in the model, which entail strong mutual dependencies between carbon and energy metabolism, protein synthesis, and growth. Figure S4 shows the mapping from resource allocation strategies to rate-yield phenotypes in the case of batch growth on glucose. 10 Whereas every strategy gives rise to a unique rate-yield phenotype (Methods), the inverse is not true: infinitely many strategies can account for a given growth rate and growth yield. Intuitively, by Eqs. 1 and 2, a given rate-yield phenotype fixes the substrate uptake rate v mc and the sum v mer + ρ mef v mef + (ρ rq -1) (v r + v q )), representing the loss of carbon due to CO 2 outflow and acetate secretion. Different protein concentrations can lead to fluxes that add up to the latter sum, and thus enable the cells to grow at the specified rate and yield (Figure S4).

How do the resource allocation strategies underlying a given rate-yield phenotype compare with those actually observed? Unfortunately, for most of the strains in the database, other than the BW25113 strain used for calibration, no proteomics data are available. Exceptions are the other commonly-used wild-type strains MG1655 and NCM3722 [START_REF] Schmidt | The quantitative and condition-dependent Escherichia coli proteome[END_REF]. Schmidt et al. found growth rates equal to 0.67 h -1 and 1.03 h -1 for MG1655 and NCM3722, respectively [START_REF] Schmidt | The quantitative and condition-dependent Escherichia coli proteome[END_REF]. Growth yields reported by other groups for these strains, in the same conditions with approximately the same growth rates, are 0.54 for MG1655 and 0.6 for NCM3722 [START_REF] Long | Fast growth phenotype of E. coli K-12 from adaptive laboratory evolution does not require intracellular flux rewiring[END_REF]48]. We compared resource allocation strategies predicting the above rates and yields with those actually observed for MG1655 and NCM3722 (Methods).

As seen in Figure 4A, there exist strategies predicting growth at approximately the same rate and yield as MG1655 that fall within a 10% margin of the experimental values. This indicates that the model accounts for the observed strategy of a strain with a similar growth phenotype as the BW25113 strain. NCM3722 grows much faster and at a higher yield than BW25113. Whereas there exist strategies with a good correspondence for the ribosomal protein fraction φ r , the predicted φ c values for a strain with the same rate and yield as NCM3722 are much higher than those observed (Figure 4B). In other words, the model predicts a higher protein fraction for enzymes in central carbon metabolism (m /p) than observed in the proteomics data. The underlying problem is that in our model the carbon uptake and metabolization rate is directly proportional to the enzyme concentration (v mc = k mc m c S/(S + K mc ), where k mc [h -1 ] is an environment-dependent apparent catalytic constant and S the constant substrate concentration, Text S1). Therefore, an increase in v mc necessary for the higher growth rate of the NCM3722 strain, requires a proportional increase in enzyme concentration, and causes the protein fraction m c /p to increase. The measured protein fraction, however, remains approximately stable (0.09 for BW25113 and 0.07 for NCM3722). A similar but opposite divergence of model and data is seen in the case of the protein fraction of enzymes in energy metabolism (Figure 4B).

Resource allocation explains the protein fraction of ribosomes necessary to attain the observed rate and yield. The discrepancies between predicted and observed enzyme fractions, however, show that additional regulatory factors need to be considered in the latter case. This conclusion agrees with the widespread view that the regulation of fluxes in central metabolism is largely due to the regulation of enzyme activity rather than enzyme concentrations [START_REF] Davidi | Lessons on enzyme kinetics from quantitative proteomics[END_REF][START_REF] Donati | Crosstalk between transcription and metabolism: how much enzyme is enough for a cell?[END_REF]. While little is known about the mechanisms allowing NCM3722 to grow much faster than BW25113, genomic changes and their physiological impact have been identified for ALE strains [START_REF] Lacroix | Use of adaptive laboratory evolution to discover key mutations enabling rapid growth of Escherichia coli K-12 MG1655 on glucose minimal medium[END_REF][START_REF] Utrilla | Global rebalancing of cellular resources by pleiotropic point mutations illustrates a multi-scale mechanism of adaptive evolution[END_REF][START_REF] Cheng | Global metabolic network reorganization by adaptive mutations allows fast growth of Escherichia coli on glycerol[END_REF]. In an ALE mutant evolved in glycerol, the change in growth rate was attributed to a change in activity of the GlpK enzyme [START_REF] Cheng | Global metabolic network reorganization by adaptive mutations allows fast growth of Escherichia coli on glycerol[END_REF], leading to higher glycerol uptake rates. In the model, the latter mutation would translate to an increase in the catalytic constant k mc (Text S1).

In order to verify the hypothesis that changes in kinetic parameters can restore the consistency between predicted and observed enzyme concentrations, we modified the analysis of 11 the model. Instead of varying only resource allocation parameters (φ q , φ r , φ c , φ er , φ ef ), we also allowed the catalytic constants (k mc , k mer , k mef ), representing the (apparent) enzyme turnover rates in central carbon and energy metabolism (Text S1), to increase or decrease by at most a factor of 2. The results of the simulations are shown in Figure 4C-D. For both the MG1655 and NCM3722 strains, there now exist resource allocation strategies and corresponding protein fractions consistent with the experimental data within a 5% margin. The NCM3722 strategies display an increased value of k mc (Figure S5). That is, the model predicts that glycolytic enzymes are more active in NCM3722 as compared to BW25113 during growth on glucose. As a consequence, resources can be shifted from glycolytic enzymes to other proteins and metabolites, which allows the growth rate and yield to increase. While resource allocation provides quantitative bounds on growth rate, growth yield, and uptake and secretion rates agreeing with the experimental data, quantitatively correct predictions of protein concentrations therefore require metabolic regulation to be taken into account as well.

Discussion

Analysis of the resource allocation strategies adopted by microbial cells can explain a number of phenomenological relations between growth rate, growth yield, and macromolecular composition [START_REF] Scott | Interdependence of cell growth and gene expression: Origins and consequences[END_REF][START_REF] Scott | Emergence of robust growth laws from optimal regulation of ribosome synthesis[END_REF]4,[START_REF] Giordano | Dynamical allocation of cellular resources as an optimal control problem: Novel insights into microbial growth strategies[END_REF][START_REF] Weiße | Mechanistic links between cellular tradeoffs, gene expression, and growth[END_REF][START_REF] Reimers | Cellular trade-offs and optimal resource allocation during cyanobacterial diurnal growth[END_REF][START_REF] Bosdriesz | How fast-growing bacteria robustly tune their ribosome concentration to approximate growth-rate maximization[END_REF][START_REF] Towbin | Optimality and suboptimality in a bacterial growth law[END_REF][START_REF] Maitra | Bacterial growth laws reflect the evolutionary importance of energy efficiency[END_REF][START_REF] Dourado | An analytical theory of balanced cellular growth[END_REF]. We have generalized this perspective to account for a striking observation: the large variability of rate-yield phenotypes across different strains of a bacterial species grown in the same environment. We constructed a coarse-grained resource allocation model (Figure 1), which was calibrated using literature data on batch and continuous growth of E. coli in minimal medium with glucose or glycerol. In each of the conditions, we considered the space of rate-yield phenotypes predicted by the model when allowing resource allocation to vary over the entire range of possible strategies, while keeping the kinetic parameters constant.

This approach is based on a number of strong assumptions. The coarse-grained nature of the model reduces microbial metabolism and protein synthesis to a few macroreactions, instead of accounting for the thousands of enzyme-catalyzed reactions involved in these processes [48,[START_REF] Adadi | Prediction of microbial growth rate versus biomass yield by a metabolic network with kinetic parameters[END_REF][START_REF] Wortel | Metabolic enzyme cost explains variable trade-offs between microbial growth rate and yield[END_REF]. Resource allocation is reduced to constraints on protein synthesis capacity, whereas other constraints such as limited solvent capacity and membrane space may also play a role [START_REF] Adadi | Prediction of microbial growth rate versus biomass yield by a metabolic network with kinetic parameters[END_REF][START_REF] Beg | Intracellular crowding defines the mode and sequence of substrate uptake by Escherichia coli and constrains its metabolic activity[END_REF][START_REF] Zhuang | Economics of membrane occupancy and respiro-fermentation[END_REF][START_REF] Szenk | Why do fast-growing bacteria enter overflow metabolism? testing the membrane real estate hypothesis[END_REF]. All possible combinations of resource allocation parameters were considered, limited only by the constraint that they must sum to 1. Observed variations in protein abundance are less drastic [START_REF] Schmidt | The quantitative and condition-dependent Escherichia coli proteome[END_REF][START_REF] Hui | Quantitative proteomic analysis reveals a simple strategy of global resource allocation in bacteria[END_REF], and coupled through shared regulatory mechanisms [START_REF] Scott | Emergence of robust growth laws from optimal regulation of ribosome synthesis[END_REF][START_REF] Chubukov | Coordination of microbial metabolism[END_REF]. The kinetic in the model have apparent values absorbing unknown regulatory effects, specific to each growth condition. This contrasts with strain-specific kinetic models with an explicit representation of the underlying regulatory mechanisms [START_REF] Weiße | Mechanistic links between cellular tradeoffs, gene expression, and growth[END_REF][START_REF] Erickson | A global resource allocation strategy governs growth transition kinetics of Escherichia coli[END_REF][START_REF] Millard | Metabolic regulation is sufficient for global and robust coordination of glucose uptake, catabolism, energy production and growth in Escherichia coli[END_REF], and does not allow our model to be used for transitions between growth conditions. Despite these limitations, we observed a very good quantitative correspondence between the predicted and the measured rate-yield phenotypes across different E. coli strains grown in the same environment (Figure 2). This suggests that differences in resource allocation are a major explanatory factor for the observed rate-yield variability. We also found a surprising absence of correlation between growth rate and relative acetate overflow (Figure 3 and Figure S1), contrary to what has been observed for a single strain in different conditions (Figure S3). This observation, A. Observed research allocation strategy for batch growth of the MG6155 strain in minimal medium with glucose [START_REF] Schmidt | The quantitative and condition-dependent Escherichia coli proteome[END_REF] (red dots) as well as examples of resource allocation strategies giving a similar growth rate and growth yield in the simulations of Figure 2B (grey dots). The predicted growth yield Y , growth rate µ, and resource allocation parameters φ q , φ r , φ c , φ er + φ ef are each allowed to vary by less than 10% from the observed values. B. Idem for the NCM3722 strain, without imposing φ c and φ er + φ ef to agree with the experimental data. The actually observed strategies do not match the predicted strategies in the case of NCM3722. C. As in panel A, but allowing the apparent catalytic constants (k mc , k mer , k mef ) to vary by a factor of 2 in comparison with the reference values determined for the BW25113 strain (Text S2). The discrepancy between the predicted and observed strategies is set to a 5% limit. D. Idem for the NCM3722 strain. When relaxing the constraint of constant kinetic parameters, the model agrees with the observed resource allocation strategies of both MG6155 and NCM3722. This requires an increase of k mc and a decrease of k mef (Figure S5). [START_REF] Forchhammer | Growth rate of polypeptide chains as a function of the cell growth rate in a mutant of Escherichia coli 15[END_REF] which was made before for a selection of ALE mutant strains [48], raises several questions. First, why is a shift of resources from respiration to fermentation, reflecting a trade-off between energy efficiency and protein cost [4, [START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF][START_REF] Pfeiffer | Cooperation and competition in the evolution of ATP-producing pathways[END_REF], not required for obtaining a higher growth rate? We inspected the strategies producing a high growth rate and little acetate overflow, and compared these with resource allocation in the reference strain used for calibration. For example, high ratehigh yield strains limit their investment in ribosomes, but increase concentrations of metabolic precursors and/or energy cofactors to maintain high protein synthesis fluxes (Figure S6). This suggests a trade-off reminiscent of the hypothesized trade-off between enzyme and metabolite concentrations in central carbon metabolism [START_REF] Fendt | Tradeoff between enzyme and metabolite efficiency maintains metabolic homeostasis upon perturbations in enzyme capacity[END_REF][START_REF] Dourado | On the optimality of the enzyme-substrate relationship in bacteria[END_REF]. Second, why do most E. coli strains in Figure 3 nevertheless display acetate overflow and grow at a suboptimal yield? The answer is not known, but probably relates to the natural environment in which E. coli has evolved, characterized by periodic anoxic conditions and synthrophic interactions with other species [START_REF] Van Elsas | Survival of Escherichia coli in the environment: fundamental and public health aspects[END_REF].

When comparing the predicted resource allocation strategies with those observed for a given rate and yield (Figure 5A-B), we found some discrepancies that cannot be solely attributed to the uncertainty in the proteomics data. In particular, the concentrations of enzymes in central carbon metabolism are consistently overestimated by the model for rates and yields exceeding those observed for the BW25113 reference strain, due to the assumption that the kinetic parameters are constant. When allowing the catalytic constants of the metabolic macroreaction to vary as well, we could restore consistency between the model and the data (Figure 5). Further investigation indicated that in strains attaining a higher growth rate and yield than the BW25113 strain, glycolytic enzymes must be more active. To our knowledge, no experimental data exist to specifically test this prediction. It is known, however, that the activity of pyruvate kinase, regulated by fructose-1,6-bisphosphate [65], increases with a higher glycolytic flux [START_REF] Kochanowski | Functioning of a metabolic flux sensor in Escherichia coli[END_REF][START_REF] Kremling | Analysis of global control of Escherichia coli carbohydrate uptake[END_REF]. In other words, a higher growth rate is accomplished at a higher activity of this glycolytic enzyme.

We conclude that the assumption that resource allocation alone can explain the observed rate-yield variability is a good first approximation that can be refined by taking into account differences in enzyme activity, either due to genetic differences between strains or to regulatory mechanisms responding to physiological changes. Predictions of the range of rates and yields achievable by a microbial species, and relating these to underlying resource allocation strategies, are of great interest for a fundamental understanding of microbial growth. In addition, by extending the model with a macroreaction for the production of a protein or a metabolite of interest [START_REF] Yegorov | Optimal control of bacterial growth for the maximization of production[END_REF], this provides rapidly exploitable guidelines for metabolic engineering and synthetic biology, by pointing at performance limits of specific strains and suggesting improvements. While instantiated for growth of E. coli, the model equations are sufficiently generic to apply to other microorganisms. The calibration of such model variants can benefit from the same hierarchical procedure as developed here, exploiting largely available proteomics and metabolomics datasets.

Methods

The resource allocation models were derived from a limited number of assumptions on the processes underlying microbial growth, as explained in Text S1. The parameters in the models were determined from literature data, as described in Text S2. In order to produce the rate-yield phenotype plots, we uniformly sampled combinations of resource allocation parameters φ r , φ c , 14 φ er , and φ ef such that their sum equals 1-φ q , where φ q was set to a constant value determined from the data (Text S2). Starting from initial conditions, the system was simulated for each combination of resource allocation parameters until a steady state was reached. When sampling the space of initial conditions for several resource allocation strategies, the system was found to always reach the same steady state. All simulations were carried out by means of Matlab R2020b. The models and the simulation code used for generating all figures in the paper are available as File S4.

The rate-yield database was compiled from the experimental literature (Files S1 and S2). When yields were explicily reported, they were converted to the dimensionless quantity Cmmol substrate Cmmol -1 biomass used in this paper, using appropriate conversion constants (Text S1). The relative acetate secretion flux was computed for those studies in the database in which both the glucose uptake rate and the acetate secretion rate were measured. Both rates were converted into the appropriate unit (Cmmol gDW -1 h -1 ).

The comparison of predicted and observed resource allocation strategies for the BW25113, MG1655 and NCM3722 strains was carried out by means of the proteomics data of Schmidt et al. ([18], Table S11). We computed the mass fraction for each protein category distinguished in the model by associating the latter with specific COG groups (r/p → Amino acid transport and metabolism and Translation; m c /p → Carbohydrate transport and metabolism; (m er + m ef )/p → Energy production and conversion; q/p → All other COG groups). The mass fraction of enzymes in energy metabolism was further subdivided into fractions attributed to respiration and fermentation, m er /p and m ef /p, in the same way as for model calibration, by distinguishing enzymes specific to fermentation, enzymes specific to respiration, and enzymes shared between respiration and fermentation ( [START_REF] Basan | Overflow metabolism in Escherichia coli results from efficient proteome allocation[END_REF] and File S3). The resource allocation strategy (φ q , φ r , φ c , φ er , φ ef ) was equated with the corresponding mass fractions.

Figure 1 :

 1 Figure 1: Coarse-grained model of microbial growth with coupled carbon and energy fluxes. Upper left figure: schematic outline of the model, showing the biomass constituents and the macroreactions, as wel as the growth and degradation of biomass. Green boxes: system of differential equations describing the carbon and energy balances, growth rate and growth yield, and resource allocation. The kinetic expressions for the reaction rates can be found in Text S1. The growth rate and growth yield are defined in terms of the fluxes of the macroreactions. Lower right figure: different protein categories considered in resource allocation. The fluxesv r , v q , v mc , v mer , v mef , v d [Cmmol or mmol gDW -1 h -1 ],the variables p, r, q, m c , m er , m ef , c, u, a * [Cmmol or mmol gDW -1 ], the resource allocation parameters φ q , φ r , φ c , φ er , φ ef [dimensionless], the degradation rate constant γ [h -1 ], the biomass density β [Cmmol gDW-1 ], the ATP yield and cost factors n mer , n mef , n r , n q [mmol Cmmol-1 ], and the CO 2 loss factors ρ mef , ρ rq [dimensionless] are formally defined in Text S1. The values of the parameters are derived in Text S2.

Figure 2 :

 2 Figure 2: Predicted rate-yield phenotype space and comparison with experimental data. A. Measured proteome fractions of the protein categories in the model (resource allocation strategies) for the BW25113 reference strain used for calibration ([18], see Text S2). B. Predicted and observed combinations of growth rate and growth yield for steady-state batch growth of E. coli in minimal medium with glucose. The rate-yield phenotypes concern the reference strain, other wild-type strains, mutant strains obtained by directed mutagenesis, and mutant strains from ALE experiments. C. Idem for batch growth of E. coli in minimal medium with glycerol. D. Idem for continuous growth in a chemostat of E. coli in minimal medium with glucose at different dilution rates (0.2, 0.35, and 0.5 h -1). The predicted yields are shown for the dilution rates ±10%. All predictions were made using the model in Figure1, calibrated for the different growth conditions, and varying the resource allocation parameters between 0 and 1, except φ q (40,000-90,000 samples). The measurements of rate and yield reported in the source literature have been converted to units h -1 (growth rate) and a dimensionless unit corresponding to Cmmol biomass Cmmol -1 substrate (growth yield) (see Text S1 and Files S1 and S2 for details).

8

 8

Figure 3 :

 3 Figure 3: Predicted flux phenotypes and comparison with experimental data. A.Glucose uptake rate (v mc ) for every predicted rate-yield combination in the case of batch growth of E. coli on minimal medium with glucose (Figure2B). B. Idem for measured combinations of glucose uptake rate, growth rate, and growth yield (File S1). The predicted and observed secretion rates are strongly correlated (R 2 = 0.93, FigureS1), in agreement with Eq. 2. C. Relative acetate secretion rate (v mef /v mc ) for every predicted rate-yield combination. D. Idem for measured combinations of glucose uptake rate, acetate secretion rate, growth rate, and growth yield (File S1). The predicted and observed relative acetate secretion rates are correlated (R 2 = 0.66, FigureS1). There is no evidence in the data of a correlation between growth rate and relative acetate secretion rate (R 2 = 0.08, FigureS1). In addition to the relative acetate overflow, the model also captures the absolute acetate overflow very well (FigureS2).

Figure 4 :

 4 Figure 4: Validation of predicted resource allocation strategies by means of proteomics data.A. Observed research allocation strategy for batch growth of the MG6155 strain in minimal medium with glucose[START_REF] Schmidt | The quantitative and condition-dependent Escherichia coli proteome[END_REF] (red dots) as well as examples of resource allocation strategies giving a similar growth rate and growth yield in the simulations of Figure2B(grey dots). The predicted growth yield Y , growth rate µ, and resource allocation parameters φ q , φ r , φ c , φ er + φ ef are each allowed to vary by less than 10% from the observed values. B. Idem for the NCM3722 strain, without imposing φ c and φ er + φ ef to agree with the experimental data. The actually observed strategies do not match the predicted strategies in the case of NCM3722. C. As in panel A, but allowing the apparent catalytic constants (k mc , k mer , k mef ) to vary by a factor of 2 in comparison with the reference values determined for the BW25113 strain (Text S2). The discrepancy between the predicted and observed strategies is set to a 5% limit. D. Idem for the NCM3722 strain. When relaxing the constraint of constant kinetic parameters, the model agrees with the observed resource allocation strategies of both MG6155 and NCM3722. This requires an increase of k mc and a decrease of k mef (FigureS5).

Acknowledgments

This work was supported by the ANR project Maximic (ANR-17-CE40-0024). The authors would like to thank Francis Mairet and Antrea Pavlou for comments on a previous version of the manuscript.

Supporting information

Text S1. Model equations (pdf file).

Text S2. Model calibration (pdf file).

Figure S1-S6. Supplementary figures (pdf file).

File S1. Database with reported rate-yield pairs for E. coli grown on glucose minimal medium (excel file).

File S2. Database with reported rate-yield pairs for E. coli grown on glycerol minimal medium (excel file).

File S3. Classification of energy proteins (excel file).

File S4. Matlab code (zip file).