isotracer: An R package for the analysis of tracer addition experiments

Tracer addition experiments, particularly using isotopic tracers, are becoming increasingly important in a variety of studies aiming at characterizing the flows of molecules or nutrients at different levels of biological organization, from the cellular and tissue levels, to the organismal and ecosystem levels. We present an approach based on Hidden Markov Models (HMM) to estimate nutrient flow parameters across a network, and its implementation in the R package isotracer. The isotracer package is capable of handling a variety of tracer study designs, including continuous tracer drips, pulse experiments, and pulse-chase experiments. It can also take into account tracer decay when radioactive isotopes are used. To illustrate its use, we present three case studies based on published data and spanning different levels of biological organization: a molecular-level study of protein synthesis and degradation in Arabidopsis thaliana, an organismal-level study of phosphorus incorporation in the eelgrass Zostera marina, and an ecosystem-level study of nitrogen dynamics in Trinidadian montane streams. With these case studies, we illustrate how isotracer can be used to estimate uptake rates, turnover rates, and total flows, as well as their uncertainty. We also show how to perform model selection to compare alternative hypotheses. We conclude by discussing isotracer’s further applications, limitations, and possible future improvements and expansions.

Tracer addition experiments are an increasingly common tool used to answer a wide variety of biological, ecological and evolutionary questions. These experiments consist 30 in injecting a labelled element into a biological system and tracing its fate throughout the system compartments across time to estimate flows of material across the com-48 ponent via a system of differential equations describing the flows of material across connected compartments, which can be outside the typical statistical expertise of biologists 50 or field geologists. While some sophisticated frameworks exist to analyze pharmacokinetics data (Gelman et al. (1996), Lunn et al. (2002) for PKBugs, Gillespie et al. (2021) 52 for Torsten, currently in development) or metabolomic data (Weindl et al., 2015), other disciplines such as evolutionary biology, ecology or ecosystem science suffer from a lack 54 of a standard analytical tool to make statistically rigorous inferences from tracer additions. For example, analysis of tracer addition experiments to study food web dynamics 56 is usually done by fitting mass balance equations focusing on one trophic compartment at a time, and trophic relationships are assumed to be known a priori, including the 58 relative proportions of diet sources when a consumer feeds on more than one source (Dodds et al., 2000;Mulholland et al., 2000;Collins et al., 2016). Such an approach 60 does not take into account the uncertainty in assumed trophic relationships (Ainsworth et al., 2010;Dodds et al., 2014), and does not allow either the uncertainty in flow esti-62 mates at one trophic level to be taken into account when estimating flow in the rest of the trophic network. Statistically reliable estimates of parameter uncertainty are par-64 ticularly important when performing comparative or experimental studies (e.g. Whiles et al. (2013); Collins et al. (2016); Norman et al. (2017); Tank et al. (2018)). 66 We developed the isotracer R package in an effort to make a statistically rigorous framework for such analyses more easily accessible for researchers. Our package system of differential equations governing the material flows, and by comparing the expected compartment sizes and tracer contents to observations. The current version 74 of isotracer implements first-order kinetics for material transfer across compartment, allows replicated units in an experimental design and can take into account categorical 76 covariates to model treatment effects. By using a Bayesian framework, we can make statistically rigorous inferences, estimate parameter uncertainty and calculate derived 78 parameters such as network-wide properties (e.g. total nutrient flow in an ecosystem). Model comparison is possible when several plausible models exist. Additionally, the 80 package also allows researchers to perform power analyses when designing their tracer addition experiments in order to make the most out of such cost-and effort-intensive 82 studies.
We first describe the mathematical framework used in our modelling approach and 84 an overview of its implementation in the package. We then illustrate the use of the package through three case studies based on published datasets: protein turnover in 86 Arabidopsis thaliana leaves (Li et al., 2017), phosphate uptake in Zostera marina individuals (McRoy and Barsdate, 1970), and nitrogen flows in a Trinidadian mountain 88 stream (Collins et al., 2016). Detailed tutorials to reproduce those case studies are available as an appendix.

90
2 Mathematical framework The system of interest can be represented as a network of connected compartments 92 ( Figure 1). As examplified in our case studies, compartments can be anything from pools of molecules (e.g. intermediates in cellular biosynthesis pathways), to tissues (e.g. 94 blood, liver, and muscle; or leaf, stem, and root), to species or functional groups (e.g. algae, invertebrate grazer, fish predator). Each compartment represents a distinct pool 96 of matter of any given chemical element of interest, and that matter can flow from one compartment to another when compartments are connected. In a tracer addition 98 experiment, a very small fraction of labelled or marked element of interest (the tracer) is added. The marked fraction is followed in time throughout the system's compartment 100 in comparison to the natural unmarked fraction. In essence, modelling data from a tracer addition experiment consists of comparing the expected trajectories of tracer 102 with observations, given a set of parameter values. In the case of stable isotope studies, the tracer is often a naturally rare heavy isotopic form of the element of interest (e.g.
x 1 x 2 x 3 y 1 λ 1 x 1 υ 3,1 x 1 υ 2,1 x 1 υ 3,2 x 2 λ 2 x 2 λ 3 x 3 values is detailed in the sections below, and can be divided into two parts. The first part is the calculation of the expected movement of material across the network for 118 those parameter values, which allows to predict the latent state of the network at any point in time. The second part is the likelihood calculation through the incorporation 120 of process and observation error: observations are assumed to be generated via some statistical distribution parameterized with the expected latent states of the network 122 calculated previously.

124
The mathematical framework described here is similar to the one described in López-Sepulcre et al. (2020). We try to keep mathematical notations consistent with it as 126 much as possible. While the presentation in López-Sepulcre et al. (2020) focuses on discrete-time modelling, we present only the equivalent continuous-time approach here, 128 which allows to specify the model based on a system of differential equations. Let's consider a network with N distinct compartments (N = 3 in Figure 1). Each 130 compartment is a pool of the element of interest (e.g. sulphur), and we can define the network state at time t by a N × 1 vector the quantity of material in compartment i at time t. Flows between connected compartments are assumed to follow first-order kinetics, and the transfer rate coefficients 134 are contained in a matrix Υ (capital upsilon) where each coefficient υ i,j determines the transfer rate from compartment j to i as υ i,j x (t) j . Transfer rates can represent a 136 variety of processes depending on the modelled system, including: chemical transformations, resource allocation among tissues, nutrient uptake, or consumption among 138 trophic levels. In addition to flowing between compartments, matter can be lost from coefficient λ i which also defines first-order kinetics for the material exiting the system from this compartment. Finally, some exogenous input can add matter to the system, 142 and is defined for each compartment i as an input function y (t) i giving the input rate for this compartment. Exogenous inputs can represent a variety of scenarios such as 144 nutrients entering an ecosystem from adjacent systems (e.g. upstream in rivers, terrestrial or aerial inputs in aquatic ecosystems), available food or nutrients to be eaten or 146 uptaken by animals or plants, or pollutants in an organism or ecosystem, to name a few examples.

148
The evolution of such a network system is entirely described by solving for x (t) in a linear system of first-order differential equations, for a given x (t 0 ) defining the initial 150 conditions. For the simple network example shown in Figure 1, the system of equations is: (1) which can be written in a matrix form: In this form, non-zero υ i,j outside the diagonal of the matrix defines a transfer from j to i, while the diagonal coefficients are the overall turnover rate coefficients for each compartment. The added vector y (t) represents the exogenous input rates.

156
In the more general case, the system of differential equations can be written: with transfer matrix A = Υ − K · I N where K is a vector containing the turnover rate 158 coefficients for each compartments such that, for compartment j, k j = λ j + N i=1 υ i,j , and I N is the identity matrix.

Likelihood calculation
As mentioned above, tracer addition data usually comprise time series of the com-162 partment sizes as well as their labelled proportions. We use interchangeably the terms "compartments" and "pools", and the terms "tracer", "labelled" element and "marked" 164 element to refer to the tracer being added, and "unlabelled" or "unmarked" element for the corresponding common isotope being measured at the same time as the tracer.

166
To model compartment sizes and labelled proportions, it is necessary to distinguish between two subpopulations of isotopes: the marked population (usually the heavy iso- giving the marked and unmarked quantities of atoms, respectively, for each compart-172 ment. These vectors are related to the total pool sizes by x (t) = m (t) + n (t) and the compartment labelled fractions are then defined as is the 174 element-wise division. In a typical tracer addition experiment, the observed time series are for x (t) and z (t) (not necessarily at the same time points), but for radioactive 176 elements the observed time series can be simply m (t) . The system of differential equations governing m (t) and n (t) is the same as for 178 x (t) above, the only differences being that y (t) is split into y m (t) and y n (t) (the vectors defining the input of marked and unmarked matter, respectively, for each compartment) 180 and that initial conditions usually differ for m (t) and n (t) . Additionally, for radioactive tracers, the λ rate coefficients are adjusted to take into account the radioactive decay 182 rate.
For a given set of parameter values, once the system of differential equations is time t, one can assume a truncated Normal distribution with the mean being the expected compartment size and a coefficient of variation cv estimated by the model: and for the observed labelled proportions one can assume, for example, a Beta distribution with mean the expected compartment labelled proportion and a precision 192 parameter φ estimated by the model: Other distributions can be a reasonable choice for the generation of observations from proportions are well below 1, and are offered as an option in isotracer.

198
The mathematical model described above can be adjusted to reflect specific features of a given tracer addition setup. i .

210
Some source compartments might be considered to be in a steady state on the time scale of the experiment. This can be the case of dissolved nutrients in a stream which 212 are constantly renewed by the water flow, or when the source compartment is extremely large relative to the amount of material which is transferred to consumer compartments 214 (e.g. ants feeding on a large source of prepared medium). To model compartment i as being in a steady state, one needs to set all the coefficients of the corresponding i th 216 row of the A matrix to zero and y i (t) to zero: this will result in a constant value for

Overenrichment
Some sampled compartments can appear over-enriched compared to their source com-220 partments. Over-enrichment happens when a receiving compartment appears more labelled than its source compartment. This occurs notably in ecosystem studies. Over-222 enrichment should not be possible according to the model, which assumes that the marked material is instantaneously mixed in the receiving compartment pool and that 224 transfer from one compartment does not preferentially affect marked or unmarked material. However, some sampled compartments might actually represent the average of 226 several sub-compartments, for example an active compartment involved in material flows and a refractory one which behaves as if isolated from the rest of the network on 228 the time scale of the experiment. This is the case, for example, for a detrital compartment sampled on a stream bed, and which might represent a mix of algae and bacteria 230 (the active portion), assimilating dissolved nutrients and preferentially grazed by invertebrate consumers, and a refractory portion of sediment and slowly decaying matter 232 (e.g. wood and leaves), not involved in significant nutrient cycling in the span of the experiment. Invertebrate consumers can appear over-enriched compared to the aver-234 age detrital compartment, but they are actually feeding mostly on the rich microbial fraction. For a split compartment i, comprised of an active and a refractory portions, 236 modelling can be adjusted by adding a π i parameter representing the active fraction of its pool at t 0 , and by only taking into account this active sub-pool when calculating 238 material flows, but by taking into account the inactive, refractory portion when calculating expected total pool sizes and apparent labelled proportions (López-Sepulcre 240 et al., 2020). 3.1 Overview The model exposed above is implemented in isotracer using a Bayesian approach.

244
The Stan program (Carpenter et al., 2017) and its R interace (the rstan package, Stan Development Team (2020) trajectories for a set of parameter values (using the Euler method with constant step size dt to numerically solve the system of differential equations), then the likelihood is 252 calculated by comparing observations to those expected trajectories. Note that, while the package user can choose the dt time step used during integration, no extensive 254 check of numerical accuracy is performed during the numerical solving, and neither the ODE solvers nor the matrix exponential tools provided by Stan are used. This 256 approach allows for greater speed of model evaluation and is expected to perform robustly when a correct dt time step and reasonable parameter priors are chosen, so 258 that only a small fraction of the material from each compartment is transferred during a time step. In practice, a good rule of thumb is to try to keep the product between dt

Interface overview
The isotracer package separates the process of modelling tracer additions into two 266 steps: (1) model definition, and (2) model fitting. This is due to the the amount of information required to define a network model, namely: network topology, initial 268 conditions, observations, and priors. Figure 2 shows an overview of the typical workflow when using isotracer and

A minimum model
The usage of the functions that define a network model (Table 1)  is applied in the model.   allowing to define versatile "pulse" designs, this function can also be used to define "drip" designs when used along set_steady(): a pulse event applied to a 336 steady-state compartment will result in a new steady state for this compartment, and "drip-on" and "drip-off" phases can thus be defined using a succession of pulse 338 events to adjust the steady state of a compartment. set_prior() to set prior distributions for individual model parameters. Implemented priors are half-Cauchy, Normal, Uniform, and scaled Beta distribution.

352
All priors are truncated to 0. Additionally, a "constant" prior is available to fix the value of some parameters during MCMC.

MCMC and posterior analysis
Once a network model is defined as explained above, MCMC sampling is performed 356 with the run_mcmc() function. This function returns an object of class mcmc.list as implemented by the coda package, and which can be used directly by many R 358 packages designed for Bayesian analyses (e.g. bayesplot). Optionally, run_mcmc() can be called with the stanfit = TRUE argument. In this case it will return the 360 raw stanfit object produced by Stan. This is especially useful when Stan produces errors or warnings during the MCMC sampling: the stanfit object can be used for diagnostics (e.g. with the shinystan package) as explained in the Stan documentation (Carpenter et al., 2017). Once a model runs without any error or warning from Stan, 364 the user is expected to use the default output from run_mcmc() (the mcmc.list object). Table 2   In addition to derived parameters and posterior predictions, isotracer can be used to calculate network properties in the particular case of a network admitting equilibrium 386 states. A network can admit an equilibrium state if all λ parameters are set to 0 (no material is lost from the network system), or if at least one compartment is set to a 388 steady state and not all λ are set to zero (a balance can be reached between network inflows and out-flows). Equilibrium states are calculated using eigenvectors of the transfer 390 matrix. The calculation is performed by the tidy_steady_states() function, which will return the compartment size estimates at equilibrium for each MCMC iteration.

392
Similarly, the tidy_flows() function can calculate the flows at equilibrium for each iteration, or the average flows over the experiment duration if the network does not 394 admit an equilibrium state. More generally, the tidy_trajectories() function will return the full compartment trajectories for each MCMC iteration if they are needed 396 for any downstream analysis.
Finally, since the package can predict compartment trajectories for a given set of pa-398 rameter values, it can also be used to simulate experiments and generate datasets. The package provides functions to sample parameter values from their priors (sample_params()), 400 to calculate trajectories for a set of parameter values (project()) and to generate predicted datasets including measurement and sampling error (sample_from()). Those 402 functions can be used to optimize experimental designs before running costly and timeconsuming real-life experiments, or to check parameter identifiability for a given design.

Case studies
The following case studies demonstrate the use of the isotracer package. In the 406 first case study on Arabidopsis protein turnover, we introduce basic concepts such as the use of covariates, and the calculation of derived parameters. In the second case 408 study on Zostera phosphate uptake, we demonstrate how to model radioactive isotopes, and how to compare derived parameters to a reference value. Finally, the study of

416
Well-regulated protein synthesis and degradation are crucial for organism development and maintenance. In a study by Li et al. (2017), 15 N labelling was used to measure 418 the degradation rate for 1228 proteins in Arabidopsis thaliana leaves from live plants.
This experiment aimed at measuring the in vivo degradation rate of proteins in young 420 plants (21-day old) in three different leaf tissues (3rd, 5th and 7th leaves). Obtaining turnover estimates for a large number of proteins allowed the authors to examine the 422 determinants of protein degradation rates, such as protein domains and involvement in protein complexes, and to identify proteins for which degradation rates differed between 424 leaves. Here, we demonstrate how to estimate turnover rates for a handful of proteins from the original study using isotracer, and how to compare those rates across leaf fraction. In addition, a similar experiment was performed in which the seeds were grown only on natural N isotope medium (mostly 14 N) and leaves were sampled at 436 similar time points (without pooling multiple individuals). Those leaves were used in chromatography and mass spectrometry analysis to estimate relative abundance 438 changes for individual proteins. The data from Li et al. (2017) is available as a Dryad repository (Li et al., 2018) 440 and the isotracer package ships several tables derived from this dataset. Here, we only focus on a handful of proteins and estimate their turnover in Arabidopsis leaves.

442
The proteins we selected are RBCL (ribulose-bisphosphate carboxylase or "rubisco", the most abundant protein in leaves in ppm), CPN60A (chaperone involved in rubisco 444 folding), PSBA (protein of the photosystem II), THI1 (thiazole biosynthetic enzyme), and PGK1 (phosphoglycerate kinase 1). Several ways exists in which the biological 446 question at hand can be answered with isotracer. For illustrative purposes, we choose protein compartments receives nitrogen "directly" from the growth medium ( Figure  4A). This topology is replicated in three instances, one per leaf tissue, and leaf identity is used as a covariate in the model so that leaf-specific parameter values are estimated. This approach assumes that no exchange of nitrogen across proteins or across leaves 452 occurs once nitrogen is incorporated into a protein. Finally, we assume that the N content in the medium is so large that it is constant compared to the protein pools, so 454 we set it to a steady state in the model ( Figure 4A). This is equivalent to the zero-order process assumed by Li et al. (2017) for protein synthesis.

456
From a modelling perspective, the dataset from Li et al. (2017) contains all the data needed for an analysis with isotracer: a time series for compartments enrichment (la-458 belled fractions) and a time series for compartments sizes (relative protein abundances). We model the network trajectories from the instant the medium was switched. We as-460 sume that the starting labelled fractions for all protein pools are the standard 15 N abundance (0.3663%); for the medium, which is set to steady state in our model, we 462 use a dummy relative abundance of 1 and a labelled fraction of 1 since almost all its N is 15 N during the experiment. Once the MCMC sampling is complete, an important 464 step is to check that the fitted model can correctly predict the original observations using a posterior predictive check ( Figure 4B). A posterior predictive check consists in 466 generating predicted data from the parameter posterior distribution, and comparing this predicted data with the actual observations used to fit the model. For example, 468 if we use a 90% probability level, we would expect the predicted trajectory envelopes to contain about 90% of the original observations. Figure 4B shows this check for the 470 labelled proportions of some compartments.
One of the aims of the original paper was to estimate the turnover rate of as many 472 proteins as possible in the three leaf types, and compare those rates across proteins and leaves. Here, we estimated the turnover rates λ for five proteins, in three leaf types 474 (since each protein pool is not transferring nitrogen to any other compartment, its turnover rate is equal to λ) ( Figure 4C). A powerful aspect of Bayesian MCMC is that 476 we can calculate posterior distributions for derived parameters by doing calculations on the primary parameters sampled during MCMC. For example, turnover times are the 478 inverse of turnover rates, and Figure 4C illustrates posteriors for turnover time which were derived from the turnover rates. the phosphate flows in a simplified network with four compartments: sediment water ("lower water"), free water ("upper water"), roots + rhizome ("lower plant"), and leaves 496 + stem ("upper plant"). To make the model more amenable, we neglect the release of phosphorus by the plant into the water and only consider four connections in the 498 network model: from upper water to leaves and stem, from lower water to roots and rhizome, and bidirectional flow between leaves and stem (upper plant) and roots and 500 rhizome (lower plant). We examine the light effect by using the light treatment as a covariate in our model, while the addition of 32 P to either the upper or the lower 502 compartment is considered as a replication within light treatment. In other words, we assume that the rate coefficients defining the phosphorus flows do not depend on 504 the compartment where 32 P was added, but that the light conditions could have an effect on the rate coefficients. The radioactive decay of 32 P is automatically taken 506 into account by isotracer by adding the appropriate decay rate to the estimated loss rate coefficient of each compartment once the half-life has been specified in the model.

508
Note that since only 32 P (not total P) was measured in these experiments, our model considers only 32 P pools and we fixed the labelled fractions of all observations to 1 510 (fully labelled). We compiled the data from McRoy and Barsdate (1970) from the original article 512 figures and tables. We used some approximations to convert the reported values to usable time series. Notably, the reported data for each time point was in cpm/mg 514 of dry weight, and we used a single mean value of each tissue total dry weight per treatment to convert those values into 32 P quantities per compartment for the purpose 516 of this illustrative case study. Given those approximations, the results reported here should be taken as an illustration of the use of isotracer, rather than as a biologically 518 precise report on eelgrass physiology. After fitting the model to determine the transfer rate parameters, we estimated 520 the average flows between compartments over the time of the experiment, in each light treatment ( Figure 5A). The Sankey plot shown in this figure is a useful way to visualize 522 flows across a network, but it should be noted that it only includes point estimates of the flows. For a statistically robust comparison between treatments, we can perform 524 calculations on the MCMC values of the primary parameters to build posteriors for derived parameters. For example, to test if the uptake rate coefficients differed between 526 light treatments, we calculated the posteriors for their ratios between light and dark conditions ( Figure 5B). Based on Figure 5B, only the coefficient for phosphate uptake 528 from upper water to the leaves and stem increased significantly in light. For the other coefficients, the 95% credible intervals were so large that they overlapped with a ratio 530 value of 1. This large uncertainty is due to the variability of the original data and the blunt assumptions that we made when we converted the original data to usable 532 time series. We expect that a dataset where total tissue dry weight were estimated at the same time as cpm/mg dry weight would allow to reduce the uncertainty a lot.

534
However, this illustrates the ability of our method to estimate uncertainty in derived parameters and to test rigorously for differences between treatments.

Nitrogen cycle in a Trinidadian mountain stream
Tracer addition experiments have been used profitably in ecology to identify the strength 538 of trophic connections and to quantify nutrient cycling in ecosystems. In a study by Collins et al. (2016), the authors quantified the nitrogen exchanges in the foodwebs 540 of mountain streams in Trinidad, from the dissolved nutrients (NH + 4 and NO -3 ) to the invertebrate consumers. They examined the effect of light conditions and of the 542 presence/absence of a fish consumer on the nitrogen dynamics in the foodweb. Data was collected by dripping 15 N-enriched ammonium into two streams in Trinidad, and 544 samples from each foodweb compartment were taken during the drip and after the drip in several transects in each stream. The transects were located at different lo-546 cations downstream of each drip. The drip phase lasted 10 days, and the post-drip phase lasted 30 days (Figure 3). Compartment enrichment was measured by deter-548 mining 15 N/ 14 N ratios by mass spectrometry, while compartment sizes were estimated in mgN/m 2 using various field sampling techniques depending on the compartment 550 type (dissolved, primary producer, invertebrate). Streams differed in their exposure to sunlight (natural versus trimmed canopy) and transects within streams differed in the 552 presence/absence of guppies (Poecilia reticulata). The dataset used in Collins et al. (2016) is shipped with isotracer as the lalaja table. In this case study, we will use the data from only one stream and a subset of the original foodweb for simplicity, but a full analysis including both streams and using the framework used in isotracer was 556 presented in López-Sepulcre et al. (2020). We show how isotracer can be used to test for the existence of a candidate trophic link by comparing models, to determine active and refractory portions of inhomogeneous compartments, and to estimate the relative proportions of inputs into a given compartment.
We ran two models using the topology presented in Figure 6A. Those models differed only by the inclusion or not of a candidate trophic link between the grazer genus 562 Psephenus and the predator genus Argia (dashed arrow in Figure 6A). Dissolved nutrients were set to a steady state given that they were renewed by the stream flow, and 564 we allowed the two primary producer compartments, the algal cover on rock subtrate (epilithon) and the fine benthic organic matter (FBOM), to be modelled as split (in-566 homogeneous) compartments given that the field sampling cannot distinguish between active and refractory portions of those compartments. We used data from three tran-568 sects from one stream (Upper LaLaja, for which the canopy cover was thinned). The transects differed in their effective addition regime due to their different distances from 570 the drip source, but otherwise were used as simple replicates in our model (i.e. they shared the same parameter values). We used DIC to compare the models: the ∆DIC 572 between the models was relatively small (2.73), and suggested that the model without the candidate trophic link fitted the data better, with a DIC-derived weight of 80% for 574 this model (López-Sepulcre et al., 2020). Posterior predictive checks were statisfactory ( Figure 6B).

576
As for the Zostera case study, we can use a Sankey plot to visualize flows between compartments ( Figure 6C). In this case, the area of a given compartment is 578 proportional to the compartment size (a quantity of nitrogen) while the thickness of the connecting ribbons is proportional to the flows (a quantity of nitrogen per time 580 unit). It results from this that the length of a given compartment is proportional to its turnover time (per time unit). The Sankey plot clearly shows the difference in flows 582 between nutrients and primary producers, compared to between primary producers and upper trophic levels. The estimated active portions in epiliton and FBOM were 0.51 584 (95%CI 0.25-0.77) and 0.05 (95%CI 0.01-0.14), respectively. Being able to estimate an active portion instead of having to set it to a fixed value a priori is important when no 586 prior knowledge is available to make a robust, educated guess about this value. It is also the case for proportions of different inputs: the modelling approach implemented 588 in isotracer estimates proportions of different inputs into a compartment from the tracer data, rather than having to rely on estimates from e.g. gut content analysis or 5 Concluding remarks: current limitations, exten-596 sibility, and future directions The isotracer package implements a Bayesian approach to model tracer addition 598 experiments and to estimate parameter uncertainty rigorously. Its primary aim is to allow researchers whose typical expertise might not cover statistical modelling of systems governed by ODE to perform rigorous and reliable statistical analyses, so that the valuable data produced by tracer addition experiments can be fully taken advantage of. The package allows for modelling of a variety of addition designs (pulse, drip, pulse-chase), can take into account both stable and decaying (e.g. radioactive) 604 tracers, and can be used for model comparison to test hypotheses relating to network structure. The toolkit provided for post-MCMC analyses is versatile and can be used to calculate network-wide properties (and their uncertainties) such as flow estimates or compartments steady state sizes when applicable.
From a technical point of view, the numerical solver currently used is reasonably fast (running 2000 MCMC iterations for the Trinidadian stream case study presented here 610 took 90 seconds on a 3.6 GHz processor with one core per chain). This performance is at the expense of automatic numerical accuracy checks and the solver must thus 612 be used with some caution. Stan diagnostics can help detect situations where this approach fails in obvious ways, but it is strongly recommended in all cases to perform 614 some extraneous MCMC runs with a few smaller time step values once the model behaves appropriately in order to check for numerical stability. Future versions could 616 implement an option to use Stan's ODE solver or matrix exponentials for increased robustness, at the cost of speed of execution, to facilitate this checking step.

618
From a modelling perspective, a future improvement could include the implementation of continuous covariates in addition to discrete ones. This would allow the study 620 of phenomena such as size-specific uptake in consumers, or temperature-specific resource allocation. Isotope fractionation could also be taken into account: the current 622 version of isotracer ignores it and assumes that the tracer and the corresponding "unmarked" material are governed by identical transfer rate coefficients. For ranges of 624 isotope enrichment usually observed in trophic networks (e.g. δ 15 N typically increasing by about 3.4 per trophic level (Cabana and Rasmussen, 1994), but with large vari-626 ations depending on the study system, e.g. Mill et al. (2007)), it is likely that this is a reasonable approximation given the very large enrichment caused by the experimental 628 addition (e.g. δ 15 N up to about 2500 in some primary consumers in Collins et al. (2016), but as low as 40 in higher trophic levels further away from the drip source).

630
However, to accurately take into account isotopic fractionation at higher trophic levels or due to other processes (Robinson, 2001), an option to model it when it is known to 632 occur could be implemented. Importantly, the current version of isotracer only models first-order reactions: 634 the quantity of transported material per unit of time is calculated as the product of a constant rate coefficient and of the source compartment size. Such first-order reac-636 tions result in simple systems of differential equations, and offer in most cases a good approximation to the modelled network dynamics on the time scale of the experiment.

638
However, in some scenarios, second-order reactions, or more complex equations, govern the material flows. In other cases, the rate coefficients are not constant over the 640 time scale of the experiment (e.g. time-dependent rate of glucose uptake from bloodstream after injection). The isotracer package can be extended with more complex 642 models to be able to model such cases, including cases where several sources must react together following stoichiometric constraints (e.g. kinetic flux profiling of known metabolic pathways, Yuan et al. (2006)). Such a task would involve extending both the Stan code (to calculate the expected compartment trajectories in the Stan model) and