Description of the continuous nature of organic matter in models of soil carbon dynamics

The understanding of soil organic matter (SOM) dynamics has considerably advanced in recent years. It was previously assumed that most SOM consisted of recalcitrant compounds, whereas the emerging view considers SOM as a range of polymers continuously processed into smaller molecules by decomposer enzymes. Mainstreaming these new insights in current models is challenging because of their ill-adapted framework. We propose the C-STABILITY model to resolve this issue. Its innovative framework combines compartmental and continuous modelling approaches to accurately reproduce SOM cycling processes. Model theoretical simulations highlight the role of SOM accessibility on its turnover. They reveal how enzyme depolymerization, decomposer community succession and carbon use efficiency are key drivers of SOM composition and quantity during decomposition and at steady-state. The mathematical structure of C-STABILITY may be tailored to different scales and is open to alternative formulations, thus offering a promising foundation for enhancing carbon cycling predictions.


Introduction
Soil organic matter (SOM) is the largest reservoir of organic carbon (C) on land (Lehmann and Kleber, 2015). Understanding and modeling the processes driving its dynamics is essential to predict SOM response to changes in climatic conditions and human land management, and its contribution to soil functions and climate mitigation.
Our view of the nature of SOM and decomposition pathways has recently been challenged by the wealth of information generated by new analytical techniques (Blankinship et al., 2018). It was previously assumed that most SOM consisted of recalcitrant litter material and humic molecules formed by the condensation of decaying substrates. The emerging view is now that SOM occurs as a range of organic compounds continuously processed into smaller molecules by microorganisms through the production of extracellular enzymes (Lehmann and Kleber, 2015). Under this new SOM cycling paradigm the local environment is seen a critical driver of SOM decomposition as it determines its physicochemical accessibility to extracellular enzymes and modulates microbial metabolism (Kemmitt et al., 2008;Schmidt et al., 2011;Schimel and Schaeffer, 2012;Dungait et al., 2012).
Metabolic activity of microbes consists of catabolism, in which extracellular enzymes depolymerize SOM, and anabolism, which is responsible for CO 2 emissions and biochemical transformation of assimilated compounds. The biosynthesized microbial metabolites are known to strongly interact with protective mineral surfaces and may be sequestered for a long time in the soil (Schimel and Schaeffer, 2012;Cotrufo et al., 2015;Liang et al., 2017Liang et al., , 2019. In the last decade new information on the functional diversity of microbial communities and on their catabolic action has become available. Genome analyses notably enable characterization of the enzyme sets encoded by each decomposer species and classification of decomposer taxa into functional guilds (Eastwood et al., 2011;Floudas et al., 2012;Riley et al., 2014;Kohler et al., 2015;Talbot et al., 2015;Nagy et al., 2016). Enzymes production patterns can now be monitored over time via proteome and secretome analyses (Zhang et al., 2016(Zhang et al., , 2019.
How have models of SOM dynamics evolved in recent years to incorporate novel knowledge and theories?
A number of new compartment models better representing microbial ecology processes have been proposed (Todd-Brown et al., 2012;Sulman et al., 2018). The most often incorporated features include decomposer physiological factors (substrate use efficiency, growth and mortality rates, etc.) and extracellular enzyme properties such as production rate, catalytic properties, etc. (Moorhead and Sinsabaugh, 2006;Allison, 2012;Kaiser et al., 2015). The microbial response to fluctuating environmental parameters such as temperature and nitrogen availability is also sometimes integrated (Treseder et al., 2012). Some models explicitly describe how the soil spatial microarchitecture governs decomposer access to substrates (Monga et al., 2014;Vogel et al., 2015;Kaiser et al., 2015) or specify the SOM protection mechanism Abramoff et al., 2018). Yet these compartment models are limited because they assume that SOM can be represented as a few discrete pools with differing turnover times (Blankinship et al., 2018). They do not describe the continuum of SOM at different stages of decay. Hence they cannot reproduce how enzyme functional diversity and decomposer ecology regulate C turnover through the progressive depolymerization of SOM until it becomes accessible for uptake and subsequent assimilation or mineralisation. Besides, the description of a large number of compounds and processes leads to model prediction uncertainty (Shi et al., 2018).
In contrast to compartment models, continuous models have received little attention from developers in recent years. Existing models represent SOM as a distribution along a quality axis and describe how organic matter (OM) quality evolves towards more recalcitrant -or stable -stages as decomposition progresses, without a detailed description of microbial and biogeochemical mechanisms (Ågren et al., 2017Bartsev and Pochekutov, 2019;Menichetti et al., 2019). At first glance, they seem particularly appealing for representing the current view of SOM as a range of heterogeneous organic compounds at different decomposition stages with few parameters. But they have three major shortcomings. First, quality is viewed as an intrinsic property of OM, generally in line with the outdated recalcitrance notion. Second, quality is a concept defined relative to a uniform decomposer group, which is at steady-state in existing models what is not realistic.
These two characteristics hinder the investigation of issues related to functional ecology and decomposer community succession. Third, the notion of quality is complex, not measurable and this abstraction is a major barrier to wide acceptance of the continuum approach.
Here we propose a new general model of SOM dynamics combining advantages and suppressing drawbacks of the current compartment and continuum quality models. This Carbon Substrate Targeted AccessiBIL-ITY (C-STABILITY) model is driven by microbial activity and combines compartmental and continuous approaches. The substrate is split into biochemical classes, separating pools accessible to enzymes from inaccessible ones (Figure 1a). Within each pool a continuous approach is implemented to describe the level of substrate polymerization (Figure 1b). Enzyme action and microbe uptake are targeted to specific chemical forms of the substrate, depending on its biochemistry and physicochemical accessibility. The model distinguishes between substrate accessibility to microbe uptake, only possible for lowly polymerized molecules, and substrate accessibility to enzymes regulated by OM spatial arrangement within the soil matrix or interactions with other soil components (Rowley et al., 2018). Over time the enzyme-accessible substrate is fragmented The polymerization axis is oriented from the lowest polymerization level on the left to the highest polymerization level on the right. The blue domain (Denz) identifies polymers accessible to enzymes. The red domain Du identifies monomers and small oligomers accessible to microbe uptake. The amount of C corresponds to the area below the curve. and depolymerized until it eventually becomes accessible to microbe uptake. Part of the absorbed substrate is assimilated by decomposers and its biochemistry is modified by decomposer anabolism. This is accounted for in the model by the redistribution of decomposer necromass into the various substrate biochemical pools upon the death of the decomposer (Figure 1). One of the major features of the model is the description of the enzymatic polymer breakage process. A substrate cleavage factor, α enz , is introduced for each enzyme family (e.g. cellulolytic, proteolytic, etc.) to indicate if it depolymerises the substrate by cleaving its end-members or if it randomly breaks any bond ( Figure 2). This general model enables accurate and parsimonious representation of the SOM nature and degradation pathways as they are currently understood.
It notably describes the continuum of organic forms generated by enzyme depolymerization, i.e. which no current SOM compartment model has achieved to date.  Here the microbe uptake rate is nil, which explains the substrate accumulation in Du. Two substrate cleavage factor α values are tested. αenz=1 is typical of the action of endocleaving enzymes, which randomly disrupt any substrate bond and generate oligomers, inducing a rapid shift towards the Du microbial uptake domain. αenz=5 is typical of exocleaving enzymes attacking the substrate end-members to release small oligomers, inducing a slower shift towards Du.
In this article, we explore C-STABILITY predictions by investigating four scenarii inspired by current keyquestions and designed on the basis of a wide range of recent publications. First, we examine the sensitivity of the model parameters in a simple case of cellulose decomposition to show how enzyme cleavage impacts the substrate uptake by decomposer. Second, we simulate lignocellulose decomposition to illustrate how substrate accessibility to enzyme is regulated. Third, we investigate how the succession on a substrate of different functional microbe communities affects its biochemistry during decomposition. Fourth, we run the model at steady-state in a soil continuously supplied with lignocellulose and subjected to microbial recycling to determine if microbial biomass is the dominant SOM source. Finally, we discuss how the innovative C-STABILITY framework offers an attractive opportunity to improve the description and prediction of C cycling in soil.

Results
2.1 Depolymerization pattern as a major driver of organic matter uptake by microbe We ran a first scenario to illustrate how substrate accessibility to microbe uptake is regulated by enzyme catalytic activity. A cellulose polymer subjected to the cellulolytic action of a set of enzymes was split into smaller fragments, which corresponded to a gradual decrease in the cellulose polymer size distribution towards short oligomers in the microbial uptake domain, D u (Figure 3a). With the chosen parameters (supplementary table), 95% of the cellulose-C was consumed after 123 days, with a mean residence time Amount of residual cellulose-C, C available for microbial uptake, cumulated CO2, living biomass-C and dead microbial residues-C. There is no microbial recycling in this simulation. Note that the sum of cellulose-C and living microbes-C at the beginning is equal to the sum of cumulated CO2 and microbial residues-C at the end; (c) Sensitivity analysis on the amount of residual cellulose-C (± 5% uniform variability). Variance fractions are explained by the carbon use efficiency e 0 mic , mortality rate m 0 mic , substrate cleavage factor αenz, enzyme activity rate τ0 and initial microbial-C to total-C ratio.
of 97 days, which was consistent with results of laboratory incubation studies, e.g. Skyba et al. (2013).
Note that during the first 25 days the enzyme action mostly generated oligomers that were too large to be accessible for microbe uptake (Figs. 3a and 3b), inducing a decrease in the living microbe biomass.
The influence of the parameters on the amount of residual cellulose varied over time ( Figure 3c). As the decomposition dynamics were driven by microbes, any modification affecting the size of their biomass had a marked impact on the residual cellulose. Thus the mortality rate m 0 mic and the initial C mic /C tot ratio had a strong influence at the beginning of the simulation, while the microbial carbon use efficiency e 0 mic was a sensitive parameter at the end of the simulation. The enzymatic substrate cleavage factor α enz and enzyme action rate τ 0 enz , which controlled the amount of substrate available to microbes, became crucial parameters once the microbe community started growing. The decrease in the α enz value induced a faster depolymerization pattern and made the substrate accessible for decomposer uptake more rapidly ( Figure 2).
α enz indeed indicated whether the enzyme preferentially catalysed substrate end-members (exo-cleaving enzyme -high α enz values) or randomly disrupted any substrate linkage (endo-cleaving enzyme -low α enz values). Exocleaving enzymes released monomers or dimers, which caused a gradual decrease in the polymer length ( Figure 2). Endo-cleaving enzymes efficiently shifted the polymer size distribution towards short lengths. Variations in the microbe uptake rate u 0 mic had no impact on the model sensitivity, i.e. with the chosen parameters, uptake was solely limited by the amount of available carbon in D u .
This simulation evidenced that C-STABILITY differentiates substrate transformation by enzyme and substrate uptake by decomposers, contrary to other models, which do not represent substrate polymerization and consider implicitly that both processes are simultaneous. This innovation enables to report the progressive modification of the decaying substrate by catabolism and how substrate availability to microbes uptake regulates the system dynamics.

Effect of temporary substrate inaccessibility to enzyme on litter decomposition kinetics
We ran a simulation of lignocellulose decomposition to illustrate how the spatial arrangement and interactions of biochemical moieties in a substrate can induce inaccessibility to enzyme. Indeed in wood, cellulose interaction with lignin makes it inaccessible to cellulases. Cellulose depolymerisation and utilisation can only start once the lignin barrier has been altered. This typically occurred during wood attack by woodrotting fungi producing lignolytic factors and cellulases, sometimes in a temporal sequence (Zhang et al., 2016(Zhang et al., , 2019. In this second scenario, lignocellulose was made of 76% cellulose, 24% lignin. A specific design of the initial cellulose distribution was made to represent its embedment with lignin, which made it inaccessible to cellulolytic enzymes ( Figure 4, see also the video in the Additional Information). Parameters were chosen to model a decomposition pattern, in general agreement with the findings of lignocellulose decomposition studies performed in microcosms (e.g. Kaffenberger and Schilling, 2015).  Lignin deconstruction was initiated from the beginning of the simulation but the amount of lignin-C did not decrease prior to 120 days ( Figure 4). Lignin depolymerization was very slow because of the combined effects of the low lignolytic rate τ 0 enz and a high substrate cleavage factor α enz (supplementary table). Alteration of the lignin framework nevertheless rapidly opened the wood cell wall structure (Goodell et al., 2017) and cellulose became completely accessible to cellulolytic enzymes after 25 days. This required step for lignolysis delayed cellulose decomposition compared to the previous simulation with a pure cellulose substrate (95% of the cellulose-C is consumed after 140 days versus 123 days in the previous case) - (Figure 4, video in the Additional Information).
Formally, mechanisms regulating substrate accessibility to enzymes were modeled in C-STABILITY with translation of the substrate chemistry distribution from a pool where the substrate was not accessible to another pool where it became accessible to its enzymes. The substrate remained unaltered as long as it was inaccessible to enzymes. This approach conceptually differs from that adopted in many compartment models where inaccessibility is reported by diminishing the substrate turnover rate compared to its free form (e.g. Parton et al. (1987); Moorhead and Sinsabaugh (2006); Stamati et al. (2013); Wieder et al. (2014)). In the scenario where cellulose was embedded in lignin, the translation was assumed to be a linear function of the lignolytic enzymes. The simulation of substrate desorption from a mineral surface or of the physical disruption of aggregates leading to the exposure of previously inaccessible OM could be performed by considering an instantaneous translation from the inaccessible to the accessible pool.

Interactions between decomposer communities regulate the persistence of decomposing substrate and its chemistry
The C-STABILITY framework is tailored specifically to reflect the succession of microbial players during litter decomposition and their biogeochemical functions. In a third scenario we implemented the model with two microbial functional communities, also called guilds, succeeding each other on a lignocellulose substrate.
The two microbial guilds had the same biochemistry, consisting of lipids, proteins and sugars (supplementary   table), but differed in their catabolic abilities. Recent studies showed that saprophitic fungi are the first to take action in wood decomposition. They were the main producers of enzymes involved in lignocellulose degradation. Bacteria or mycorrhizal fungi peak thereafter, breaking down and utilizing dead fungal biomass rather than plant polymers (Tláskal et al., 2016;Vorísková and Baldrian, 2013;Schneider et al., 2012;Lindahl and Tunlid, 2015). We referred to them as plant degraders and microbial residue degraders, respectively.
In our simulation, the two communities also differed in their fitness. Microbial residue degraders were more competitive than plant degraders because of their higher carbon use efficiency and lower mortality rate (supplementary table). This succession is illustrated Figure 5a.  Figure 5: Effect of the succession of two decomposer communities on the residual substrate amount and biochemistry (scenario 3). The two communities have the same biochemical signature (50% sugar, 30% lipid and 20% protein). Plant decomposers (red) produce cellulolytic and lignolytic enzymes. Microbial residues decomposers (Mic. decomposer -blue) produce enzymes that depolymerize lipids, proteins and microbial sugars. The two communities differ in their fitness. Microbial residue decomposers are more competitive than plant decomposers because of their higher carbon use efficiency and lower mortality rate (supplementary table). The effect of cheating behavior was tested as follows: (a) substrate uptake is only possible for the enzyme producer or (b) uptake is also possible for the cheater but at a lower rate than that of the enzyme producer.
first 260 days. The microbial residue degrader biomass slightly declined during this first stage. It started growing once the microbial residues were released upon the death of plant degraders, after around 200 days.
Despite the strategies implemented by decomposers to achieve privileged access to the substrate they depolymerize, e.g. antibiotics or secondary metabolites secretion (Hiscox and Boddy, 2017), Schneider et al. (2012) suggested that bacteria can also proliferate on simple carbohydrates provided by fungal enzymes. This has been referred to as cheating behavior. We simulated this behavior by changing the substrate uptake rate of the community that was not producing enzymes from zero to a positive value (supplementary table  and figure). This change had direct impacts on the substrate residue chemistry at the end of the simulation.
Indeed, when microbial residue degraders had an opportunity to opt for cheating behavior and benefit from the plant compounds fragmented by plant degraders, they developed faster than the plant degraders and rapidly dominated them (Figure 5b). The decline in fungal communities involved in plant depolymerization downregulated cellulose and lignin depolymerization and disappearance rates. Lignin-C therefore persisted at the end of the simulation (supplementary figure). Interactions between different functional decomposer communities thus seemed to be a strong driver of the residual substrate chemistry.

Contributions of plant and microbe residues to SOM at steady-state
Prediction of the contribution of microbial materials to the total carbon stock is of great interest for assessing the C sequestration potential of a system. Indeed microbial residues are currently considered as being the main precursors of stable OM (Cotrufo et al., 2015). For this purpose, we resolved the analytic formulation of the C stock and chemistry at steady-state (Eqs. (13) to (15)) while considering one microbial community, microbial residue recycling, continuous plant input and the substrate accessibility to enzymes (scenario 4).
A steady-state C stock of 1.553 g C .g −1 soil was computed with the parameters given in supplementary table. The living microbial biomass accounted for 8% of the total carbon stock (Eq. (13)), which was at the upper limit but within the range observed in soil across biomes (Fierer et al., 2009). Cellulose and lignin were the most abundant chemical substances, each representing 33% of the C stock at steady state (Figure 6a).
Microbial residues contributed to the remaining 33% as follows: 16% microbial sugar, 10% lipids and 7% proteins. Within each biochemical class, the substrate was a continuum of forms at different decomposition stages. In Figure 6a, the short peak on the right of the polymerization distribution of SOM corresponded to polymers initially entering the system as plant litter or microbial residues. OM altered by enzyme action was the main contributor of steady-state OM, as shown by the dominance of polymers with a lower degree of polymerization. There was no or almost no substrate in the D u domains due to active microbial uptake of the smallest substrate fragments. The simulation closely reflected the SOM chemistry on the forest floor or in boreal soils observed in the field (Balaria and Johnson, 2013;Kohl et al., 2018). To predict the SOM chemistry in mineral soil horizons, C-STABILITY should account for the reduction in substrate accessibility due to interactions with the mineral phases, particularly for microbial products that are frequently involved in aggregates or bind on mineral surfaces (Miltner et al., 2012;Cotrufo et al., 2015;Clemmensen et al., 2013;Kallenbach et al., 2016;Wang et al., 2017).
The sensitivity analysis we carried out for the default parameters, in a system where substrates were showing the change in total C stock and its distribution among biochemical classes with ± 50% changes in the parameters (variations in uptake rate u 0 mic , carbon use efficiency e 0 mic and mortality rate m 0 mic at row 1; variations in enzyme action rate τ 0 enz and cleavage factor αenz illustrated for lipid depolymerases, cellulolytic and lignolytic enzymes at rows 2 and 3, respectively). Note that the scale of the y-axis varies.
accessible to their enzymes, confirmed the importance of microbial processes in controlling the amount and chemistry of SOM at steady-state ( Figure 6b). As already reported (e.g. Liang et al., 2017), the carbon use efficiency parameter, e 0 mic , which regulated the balance between CO 2 production and metabolite production, was a major driver. An increase in e 0 mic raised the microbial biomass, which in turn accelerated the C turnover rate through increased depolymerization and respiration rates (Eqs. (3), (11) and (12)). This effect explained the decline in plant biochemical classes with increasing e 0 mic . The biochemical classes of the microbial signature were not impacted by e 0 mic . The accelerated SOM cycling induced by raising e 0 mic was indeed counterbalanced by a higher rate of microbial necromass input into the SOM pool (Eqs. (14) and (15)).
The C stock also exhibited a strong positive linear dependence on the microbe mortality rate m 0 mic , without any distinction among the various substrate biochemistries. In contrast, the C stock and distribution seemed to be quite insensitive to changes in the u 0 mic uptake rate for the range of tested values ( Figure 6): C uptake was indeed limited here by the amount of accessible C. Consequently, the enzyme parameters controlling the availability of substrate for microbe uptake had a marked impact on the C stock and biochemistry ( Figure 6b). An increase in the enzyme action rate τ 0 enz decreased the quantity of its substrate at steadystate. In contrast, an increase in the α enz value increased the quantity of its substrate at steady-state.
It corresponded to enhanced substrate exocleavage, which promoted the release of small fragments and preserved large polymers. These results contradicted the general belief that enzyme traits have a limited effect on SOM chemistry and quantity (Schimel and Schaeffer, 2012;Liang et al., 2017) and calls for closer consideration of the impact of catabolic processes on SOM cycling.

Discussion
C-STABILITY combines compartmental and continuous modelling approaches. This enables the representation of OM as a continuum of forms and the key processes governing its cycling, such as substrate accessibility and selective depolymerization, with a parsimonious number of parameters, which no models to date have been able to achieve. C-STABILITY breaks the abstract notion of quality of previous continuous models of SOM dynamics, which is closely associated with the recalcitrance concept and inadequate to account for several decomposers, and replaces it by an operational description of the OM chemistry in terms of biochemistry, polymerization and physicochemical inaccessibility. This major change allows to recognize that substrate access to enzymes and microbes strongly regulates OM turnover. Degradation is indeed not determined by any intrinsic molecular recalcitrance or specific decay rate as in many models. It is determined by the local environment properties, the substrate biochemistry and polymerization, and the functional diversity of microbes producing enzymes. Jointly, they regulate substrate accessibility to enzyme and its selective depolymerization prior uptake.
The two modelling philosophies inspiring C-STABILITY are usually opposed and have their own limitations for long term predictions. Weakness of the continuous approach is typical of the statistical modelling.
Despite the strong predictive ability of statistical models for stable systems, the lack of explicit description of functional processes lower the confidence in predictions for perturbed systems. In opposite, compartment approach is descriptive and process-based. But, the more the system is complex, the more process-based approach is facing an explosion of its parameter number and uncertainty, making long term prediction hazardous. In C-STABILITY, SOM compartments and C fluxes are process-based inspired while the depolymerization process is depicted statistically as a continuous process. A complete process-based approach of OM depolymerization would be hazardous due to the multiplicity of chemical reactions to take into account and the scarcity of available data. Instead, the continuous distribution reports a wave of carbon traveling through the depolymerization process ( Figure 2). The moment a polymer becomes accessible to enzymes, the carbon it contained begins to flow through less and less polymerized forms. This journey is parsimoniously described by the α enz parameter, which represents the substrate cleavage type (see Lebaz et al. (2015) for a similar approach). It is open to the introduction of competition relationships between decomposers and additional limiting factors such as water and nutrient availability. Besides, the equation describing C fluxes, such as the enzyme catalytic rate and microbe mortality, were chosen here to minimize the number of parameters. Yet other formulations could also be adopted. For example, the mortality rate formulation, which assumes linear dependence on the microbial biomass as in many microbial models, exhibits a decadal oscillatory behavior in response to C input perturbations, while demonstrating little or no sensitivity of steady-state C to changes in C input (Wang et al., 2014). To avoid this behavior, Georgiou et al. (2017) proposes a density-dependent formulation of mortality, which could be implemented in C-STABILITY. The enzyme catalytic rate could also be modulated by a saturating function such as the Michealis-Mentens equation (Schimel and Weintraub, 2003), even though this was not required here to ensure model stability, as in earlier models.

C-STABILITY does not account for the immense variety of enzymes involved in SOM decomposition
Given the urgent societal demand to improve the prediction of carbon dynamics at the global scale, one great challenge in designing microbe models is to be able to effectively achieve micro-to global-upscaling (Wieder et al., 2015). It would be pointless to sum up enzyme and microbe activity predicted in an infinite number of heterogeneous environments considering the enormous computational resources that would be required. Alternatively, we postulate that the operational description of OM chemistry (biochemistry, physicochemical accessibility, polymerization) and functional decomposer guilds could be a means to effectively bridge the spatial scales. Such information is available from the subnanoscale to the plot scale through a range of techniques (spectroscopic techniques, chemical extractions, proteomics, metabar coding, etc.).
This representation of OM could enable identification of key fine-scale processes and environmental factors, which emerge when describing OM dynamics at broader scales and are relevant for inclusion in Earth system models to improve their prediction accuracy.

Substrate description
The description of OM in C-STABILITY consists of several subdivisions. First, OM composing living organisms is separated from non-living organic substrate. Several living decomposers can be considered simultaneously (e.g. bacteria, fungi, etc.). Second, non-living OM from substrate is split into biochemical classes (cellulose or plant sugar, lignin, lipid, protein and microbial sugar in this study). Third, for each biochemical class, OM accessible to enzymatic activity is separated from OM which is inaccessible due to specific physicochemical conditions (e.g. interaction between different molecules, inclusion in aggregates, sorption on mineral surfaces, etc.). For each of these pools, a continuous description of the degree of OM polymerization (p) is provided, as a distribution (Figure 1). For each biochemical class * ( * = cellulose, lignin, lipid, etc.), the range of the possible polymerization degree is identical for both inaccessible and accessible pools. The total amount of C (g C ) in each pool is where p min * and p max * are the maximum and minimum degrees of polymerization and χ * (g C .p −1 ) is the polymerization distribution. To differentiate substrate that is inaccessible from substrate accessible to enzymes, variables χ and C are denoted χ in * and C in * , and χ ac * and C ac * , respectively. The polymerization axis is oriented from the lowest to the highest degree of polymerization (Figure 1b). A right-sided distribution corresponds to a highly polymerized substrate whereas a left-sided distribution corresponds to monomer or small oligomer forms.
Accessibility to microbe uptake is described by the interval (also called domain) D u , which corresponds to small OM compounds, monomers, dimers or trimers smaller than 600 Daltons (Lehmann and Kleber, 2015), that microbes are able to take up (in red in Figure 1b). Besides, accessibility to enzymes occurs in the D enz domain (in blue in Figure 1b). Over time the substrate accessible to enzymes is depolymerized and its distribution shifts towards D u and it eventually becomes accessible to microbe uptake.

Substrate dynamics
Three processes drive OM dynamics: (i) enzymatic activity, (ii) microbial biosynthesis, and (iii) changes in local physicochemical conditions. First, enzymes have a depolymerization role which enables transformation of highly polymerized substrate into fragments accessible to microbes. Second, microbial uptake of substrate is only possible for molecules having a very small degree of polymerization. When C is taken up, a fraction is respired and the remaining is metabolized, resulting in the biosynthesis of microbial molecules that return to the substrate upon microbe death. Each microbe has a specific signature that describes its composition in terms of biochemistry and polymerization. Third, changes in local OM conditions drive exchanges between OM accessible and non-accessible to enzymes (e.g. aggregate formation/break). All of these processes are considered on a daily basis (d).

Enzymatic activity
Enzymes are specific to biochemical classes. They are not individually reported, but rather as a family of enzymes contributing to the depolymerization of a biochemical substrate (e.g. combined action of endoglucanase, exoglucanase, betaglucosidase, etc., on cellulose will be reported as cellulolytic action). The overall functioning of each enzyme family enz is described by two parameters: a depolymerization rate τ 0 enz providing the number of broken bonds per time unit and a factor accounting for the type of substrate cleavage α enz . The term F act enz (g C .p −1 .d −1 ) represents the change in polymerization of χ ac * due to enzyme activity for all p ∈ D enz , The depolymerization rate, τ enz (d −1 ), is expressed as a linear function of microbial C biomass C mic (g C ), where τ 0 enz (g −1 C .d −1 ) is the action rate of a given enzyme per amount of microbial C. If several microbial guilds are associated with the same enzyme family, we replace the C mic term by a weighted sum of the C mass of all microbial guilds involved in Eq. (3). The K enz (p −1 ) kernel provides the polymerization change from p to p, where 1 p≤p equals 1 if p ≤ p and 0 otherwise. The α enz cleavage factor denotes the enzyme efficiency to generate a large amount of small fragments and to shift the substrate distribution towards the microbe uptake domain D u (Figure 2). α enz = 1 is typical of the action of endocleaving enzymes, which randomly disrupts any bond of its polymeric substrate and generates oligomers. The shift towards D u is slower if α enz increases.
This is characteristic of exocleaving enzymes, which attack the end-members of their polymeric substrate, generate small fragments and preserve highly polymerized compounds. To satisfy the mass balance, the kernel verifies K enz (p, p )dp = 1. Then F act enz does not change the total C mass but only the polymerization distribution (i.e. F act enz (χ, p, t)dp = 0).

Microbial biosynthesis
Each microbial group (denoted mic) produces new organic compounds with assimilated C. After death, the composition of the necromass returning to each biochemical pool * of SOM is assumed to be constant, accessible and is depicted with a set of distributions s mic, * , named signature. Each distribution s mic, * (p −1 ) describes the polymerization of the dead microbial compounds returning to the pool * . The signature is normalized to ensure mass conservation, i.e. if we note that S mic, * = p max * p min * s mic, * (p)dp, then we have * S mic, * = 1. For each accessible pool of substrate, the term F upt mic, * describes how the microbes utilize the substrate available in the microbial uptake D u domain (Figure 1b). For all p ∈ D u , where u 0 mic, * (g −1 C .d −1 ) is the uptake rate per amount of microbe C. The substrate uptake rate linearly depends on the microbial C quantity.
Depending on a C use efficiency parameter e 0 mic, * (ratio between microbe assimilated C and taken up C), taken up C is respired or assimilated and biotransformed into microbial metabolites. This induces a change in the OM biochemistry and polymerization ( Figure 1a).
Finally microbial necromass returns to the substrate pools with a specific mortality, which linearly depends on the microbial C quantity, where m 0 mic (g −1 C .d −1 ) is the mortality rate per amount of microbe C (Figure 1a).

Change in local physicochemical conditions
The distribution of a substrate inaccessible to its enzymes remains unchanged over time. A specific event changing the accessibility to enzymes (e.g. aggregate disruption or desorption from mineral surfaces) is modeled with a flux from the inaccessible to the accessible pool. Transfer between these pools is described by the F loc ac, * term for each biochemistry * , F loc ac, * (p, t) = τ ac tr χ in * (p, t) where τ ac tr, * (d −1 ) is rate of local condition change toward accessibility.
Transfer in the opposite way (e.g. aggregate formation, association with mineral surfaces) is described by the F loc in, * term, where τ in tr, * (d −1 ) is rate of local condition change toward inaccessibility.

General equations
The distribution dynamics for biochemical class * is given by, −τ enz (t)χ ac * (p, t) + Denz K enz (p, p )τ enz (t)χ ac * (p , t)dp where 1 Du (p) equals 1 if p ∈ D u and 0 otherwise. The dynamics of microbial C mic is obtained by, and the CO 2 flux (g C .d −1 ) produced by the microbes is given by, mic, * Du χ ac * (p, t)dp.

Steady-state equations
The formulation of the analytic expression for the steady-state requires several assumptions. First, we only consider accessible pools for all biochemical classes and one unique group of microbes. Second, we consider that C use efficiency and uptake parameters are identical for all biochemical classes. Third, we add a constant input term to equation 9. The input distribution is denoted ι * (g C .p −1 .d −1 ) and the total carbon input flux At steady-state, the amount of microbial carbon is, For each biochemical class * , the steady-state distribution of C is obtained by combining two terms defined by the microbial uptake domain

Scenario 1: Cellulose decomposition kinetics and model sensitivity
A first simulation was run to depict cellulose depolymerisation and uptake by a decomposer community over one year (see parameters in supplementary table). We evaluated the model sensitivity to parameters using Sobol's sensitivity analysis (Sobol, 1993) (see Section 4.4). This was performed for a small parameter variation (± 5% uniform variability) using 12000 model simulations.

Scenario 2: Effect of substrate inaccessibility to enzyme on litter decomposition kinetics
A simulation of lignocellulose (76% cellulose, 24% lignin) degradation was performed by taking into account peroxidases, which deconstruct the lignin polymer, and cellulases, which hydrolyse cellulose.
The cellulose was initially embedded in lignin and inaccessible to cellulase. Progressive transfer to the accessible pool was linearly related to the activity of lignolytic enzymes in Eq. (7), where D lig. is the domain of lignolytic activity and where the τ ac,0 tr,cell. coefficient is set at 13 g −1 C for the current illustration.
The simulation was performed over one year. Enzymatic and microbial parameters given in supplementary table were chosen to be closely in line with the litter decomposition and enzyme action observation of Kaffenberger and Schilling (2015), Lashermes et al. (2016) and Zhang et al. (2016).

Scenario 3: Effect of community succession on C fluxes and substrate biochemistry
We simulated the succession of two microbial functional communities, or guilds, on the previous lignocellulose, considering microbial residue recycling. The parameters (supplementary table) were chosen according to the microbial community succession observations of Snajdr et al. (2011), Schneider et al. (2012) and Tláskal et al. (2016. The first microbial community was specialized in plant substrate degradation, the second was specialized in the degradation of microbial residues. We referred to them as plant degraders and microbial residue degraders. Microbial residue degraders were more competitive than plant degraders because of their higher carbon use efficiency and lower mortality rate (supplementary table). Both communities had the same biochemical signature, i.e. 50% polysaccharides, 30% lipids and 20% proteins. We tested the impact of cheating as follows. Either uptake was impossible, i.e. u 0 equaled 0 for the community not involved in enzyme production, or uptake was possible but at a lower rate than the enzyme producers because the substrate fragments were released in the vicinity of the enzyme producers (supplementary table).

Scenario 4: SOM composition at steady state.
We resolved the analytic formulation of the C stock and chemistry at steady-state by considering one microbial community, continuous plant input at a rate of 4 mg C .g soil −1 day −1 (Wang et al., 2012) and microbial recycling (Eqs. (13) to (15)). To be able to specifically calculate the steady-state, cellulose substrate was not embedded in lignin but directly accessible to cellulolytic enzymes ( Figure 6). Compared to previous scenarii related to laboratory experiments, enzymatic activities were divided by five to account for the impact of in-situ climatic conditions (Balaria and Johnson, 2013). Other enzymatic and microbial parameters are given in supplementary table. To explore how the catabolic traits of enzymes and anabolic traits of microbes affect the SOM composition, we modified individual parameters by 50% and documented their effects on the steady-state amount of C and its biochemical signature.

Sensitivity analysis
Sensitivity analysis encompasses a wide range of methods aimed at determining how the uncertainty of model parameters influences the model outputs. These methods improve the understanding of the model and its design.
We considered a specific method for calculating sensitivity indices defined by Sobol (1993). The method was designed to decompose the variance of a model output, denoted y, according to the various degrees of interaction between the n uncertain parameters (x i ) i∈{1,n} . Formally, by assuming that the parameter uncertainties are independent, the model output could be expressed as a sum of functions that take parameter interactions into account, .. + f 1,...,n (x 1 , ..., x n ).
This formulation leads to the definition of several sensitivity indices. The first-order Sobol's index of each parameter is, and higher order indices are defined by S i,j = Var(f i,j (x i , x j ))/Var(y), and so on. These indices are unique, with a value of 0 to 1 and their sum equals 1.
Here we focused on the first-order indices, which are used to an increasing extent in ecological modeling (Cariboni et al., 2007;Sainte-Marie and Cournède, 2019) as they are usually sufficient to give a straightforward interpretation of the actual influence of different parameters. Sobol's indices were estimated using a Monte Carlo estimator of the variance (Cournède et al., 2013).

Implementation
The model was implemented in the Julia c language (Bezanson et al., 2012). An explicit finite difference scheme approximates the solutions of integro-differential equations with a ∆t = 0.1d. time step and a ∆p = 0.01 polymerization step. Differential equations were solved with a Runge-Kutta method.