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

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ isotracer: An R package for the analysis of tracer addition experiments © 2022 The Authors. Methods in Ecology and Evolution published by John Wiley & Sons Ltd on behalf of British Ecological Society Accepted version (Final draft) Bruneaux, Matthieu; López‐Sepulcre, Andrés Bruneaux, M., & López‐Sepulcre, A. (2022). isotracer: An R package for the analysis of tracer addition experiments. Methods in Ecology and Evolution, 13(5), 1119-1134. https://doi.org/10.1111/2041-210x.13822 2022 Methods Ecol Evol. 2022;13:1119–1134. | 1119wileyonlinelibrary.com/journal/mee3 Received: 23 August 2021 | Accepted: 3 December 2021 DOI: 10.1111/2041-210X.13822 RESEARCH ARTICLE isotracer: An R package for the analysis of tracer addition experiments Matthieu Bruneaux1 | Andrés López- Sepulcre1,2,3 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2022 The Authors. Methods in Ecology and Evolution published by John Wiley & Sons Ltd on behalf of British Ecological Society. 1Department of Biological and Environmental Sciences, University of Jyväskylä, Jyväskylä, Finland 2Department of Biology, Washington University in St. Louis, St. Louis, MO, USA 3CNRS UMR 7618, Institute of Ecology and Environmental Sciences of Paris (iEES), Sorbonne Université, Paris, France Correspondence Matthieu Bruneaux Email: [email protected] Funding information Academy of Finland, Grant/Award Number: 295941 Handling Editor: Jessica Royles Abstract 1. 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. However, performing rigorous statistical analyses to gain reliable quantitative insights from these experiments often remains challenging. 2. We present an approach based on Hidden Markov Models (HMMs) 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 pulsechase experiments. It can also take into account tracer decay when radioactive isotopes are used. 3. To illustrate its use, we present three case studies based on published data and spanning different levels of biological organization: a molecularlevel study of protein synthesis and degradation in Arabidopsis thaliana, an organismallevel study of phosphorus incorporation in the eelgrass Zostera marina and an ecosystemlevel 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. 4. The isotracer package allows researchers from a broad range of disciplines to fully take advantage of their datasets through rigorous statistical analyses. We conclude by discussing isotracer’s further applications, limitations and possible future improvements and expansions. KEYWORDS molecular labelling, nutrient flows, resource allocation, web dynamics 1120 | Methods in Ecology and Evoluon BRUNEAUX ANd LóPEZ- SEPULCRE 1 | INTRODUCTION Tracer addition experiments are an increasingly common tool used to answer a wide variety of biological, ecological and evolutionary questions. These experiments consist 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 compartments. Their applications encompass all levels of biological organization, from cells and tissues (Allen & Young, 2020; Yuan et al., 2006) to organisms (Kim et al., 2016; Williams et al., 2016), communities (Freeman et al., 2013; Shik et al., 2018), and ecosystems (Collins et al., 2016; Riis et al., 2012). Tracer additions have a long history. They have been used in seminal experiments revealing the nature of the genetic material (Hershey & Chase, 1952), and more recently in highthroughput metabolomic techniques (e.g. Li et al., 2017; Yuan et al., 2006). They also have a long history in ecosystem science (Crossley & Reichle, 1969) and are becoming a common tool to study the dynamics of nutrients or elements in a variety of ecosystems including: water dynamics in a savanna (Kulmatiski et al., 2010), carbon in coral reefs (de Goeij et al., 2013), nitrogen in deciduous forest soil (Goodale et al., 2015), the effects of reindeer urine in tundra ecosystems (Barthelemy et al., 2018) and nutrient cycling in a variety of aquatic environments (Sánchez- Carrillo & Álvarez Cobelas, 2017). However, gaining reliable quantitative insights and making robust inferences from such data requires appropriate statistical techniques, which can sometimes be challenging to develop and use. Analysis of tracer addition data typically requires models incorporating a time component via a system of differential equations describing the flows of material across connected compartments, which can be outside the typical statistical expertise of biologists or field geologists. While some sophisticated frameworks exist to analyse pharmacokinetic data (Gelman et al. (1996), Lunn et al. (2002) for PKBugs, Gillespie et al. (2021) 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 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 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 relative proportions of diet sources when a consumer feeds on more than one source (Collins et al., 2016; Dodds et al., 2000; Mulholland et al., 2000). Such an approach 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 estimates at one trophic level to be taken into account when estimating flows in the rest of the trophic network. Statistically reliable estimates of parameter uncertainty are particularly important when performing comparative or experimental studies (e.g. Barneche et al., 2021; Collins et al., 2016; Norman et al., 2017; Tank et al., 2018; Whiles et al., 2013). We created the isotracer R package in an effort to develop a statistical framework for such analyses more easily accessible for researchers. Our package implements the mathematical framework recently described in López- Sepulcre et al. (2020), which is based on a similar logic as the approach used by Wollheim et al. (1999). It considers a whole network system simultaneously and uses a Bayesian approach to incorporate experimental uncertainty into model parameter estimates. For a given set of parameter values, likelihood is calculated by solving the system of differential equations governing the material flows, and by comparing the expected compartment sizes and tracer contents to observations. The current version of isotracer implements firstorder kinetics for material transfer across compartments, allows replicated units in an experimental design and can take into account categorical covariates to model treatment effects. By using a statistical framework, we can make rigorous inferences and estimate parameter uncertainty. We choose a Bayesian approach to leverage its ability to easily perform marginalization and propagate errors from estimated parameters to derived parameters of interest such as networkwide properties (e.g. total nutrient flow in an ecosystem). Model comparison is possible when several plausible models exist. Additionally, the 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 effortintensive studies. We first describe the mathematical framework used in our modelling approach and 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 Arabidopsis thaliana leaves (Li et al., 2017), phosphate uptake in Zostera marina individuals (McRoy & Barsdate, 1970) and nitrogen flows in a Trinidadian mountain stream (Collins et al., 2016). Detailed tutorials to reproduce those case studies are available as an appendix. 2 | MATHEMATICAL FRAMEWORK The system of interest can be represented as a network of connected compartments (Figure 1). As exemplified in our case studies, compartments can be anything from pools of molecules (e.g. intermediates in cellular biosynthesis pathways), to tissues (e.g. blood, liver and muscle; or leaf, stem and root), to species or functional groups (e.g. algae, invertebrate grazer and fish predator). Each compartment represents a distinct pool 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 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 compartments in comparison to the natural unmarked fraction. In essence, modelling data from a tracer addition experiment consists in comparing the expected trajectories of tracer 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. 2H, 13C, 15N, 18O, 34S) which is used to enrich the | 1121 Methods in Ecology and Evoluon BRUNEAUX ANd LóPEZ- SEPULCRE pools of the corresponding naturally abundant unmarked form (1H, 12C, 14N, 16O, 32S). Observations usually comprise both the total pool size for each compartment (i.e. the sum of the quantities of heavy and light isotopes) and the proportion of tracer in each pool (i.e. the proportion of heavy isotope to the total amount of heavy and light isotopes, often measured as δ enrichment values; Fry, 2006). In the case where only the tracer is being tracked (e.g. a radioactive tracer such as 32P), the general framework presented above can be used by considering that the observed sizes represent the tracer alone, and that the observed marked proportions are always one. Radioactive decay is taken into account in the turnover rate of the tracer. It is important to note, however, that the method works for any traceable marker (e.g. immunolabelled molecules), not just rare isotopes. The mathematical model used to calculate the likelihood of a set of parameter 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 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 of process and observation error: observations are assumed to be generated via some statistical distribution parameterized with the expected latent states of the network calculated previously. 2.1 | General formulation of transfer equations 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 much as possible. While the presentation in López- Sepulcre et al. (2020) focuses on discretetime modelling, we present only the equivalent continuoustime approach here, which allows to specify the model based on a system of differential equations. Let us consider a network with N distinct compartments ( N=3 in Figure 1). Each 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 x (t)= ( x(t) 1,x(t) 2,…,x(t) N ) where x (t) i indicates the quantity of material in compartment i at time t. Flows between connected compartments are assumed to follow firstorder kinetics, and the transfer rate coefficients (later referred to as ‘rates’ for simplicity) are contained in a matrix Υ (capital upsilon) where each rate 𝜐i,j determines the flux from compartment j to i as 𝜐 i,jx (t) j . Transfer rates can represent a variety of processes depending on the modelled system, including: chemical transformations, resource allocation among tissues, nutrient uptake or consumption among trophic levels. In addition to flowing between compartments, matter can be lost from the system (e.g. through excretion or emigration): each compartment i has a loss rate 𝜆i which also defines firstorder kinetics for the material exiting the system from this compartment. Finally, some exogenous input can add matter to the system, and is defined for each compartment i as an input function y(t) i giving the input flux for this compartment. Exogenous inputs can represent a variety of scenarios such as 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 uptaken by animals or plants, or pollutants in an organism or ecosystem, to name a few examples. The evolution of such a network system is entirely described by solving for x(t) in a linear system of firstorder differential equations, for a given x( t0 ) defining the initial conditions. For the simple network example shown in Figure 1, the system of equations is: which can be written in a matrix form: In this form, nonzero 𝜐i,j outside the diagonal of the matrix defines a transfer from j to i, while the diagonal coefficients are the overall turnover rates for each compartment. The added vector y(t) represents the exogenous input fluxes. In the more general case, the system of differential equations can be written: (1) ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ dx1 dt=− �𝜐2,1 +𝜐3,1 +𝜆1�x(t) 1+y(t ) 1 dx2 dt=𝜐2,1x(t) 1−�𝜐3,2 +𝜆2�x(t) 2 dx3 dt=𝜐3,1x(t) 1+𝜐3,2x(t) 2−𝜆3x(t) 3 (2) ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ dx1 dt dx2 dt dx3 dt ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ −�𝜐2,1 +𝜐3,1 +𝜆1�00 𝜐2,1 −�𝜐3,2 +𝜆2�0 𝜐3,1 𝜐3,2 −𝜆3 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ x(t) 1 x(t) 2 x(t) 3 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ + ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ y(t) 1 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . (3) dx dt =Ax(t)+y(t ) FIGURE 1 Example of a basic network with three compartments. The xi values in each compartment are quantities of matter. The values along arrows are fluxes (quantity per time unit). y1 is a flux (quantity per time unit); 𝜐i,j and 𝜆i are rate coefficients of firstorder kinetics (per time unit), simply called ‘rates’ in this study for convenience. xi and y1 are functions of time t while 𝜐i,j and 𝜆i are constant rates 1122 | Methods in Ecology and Evoluon BRUNEAUX ANd LóPEZ- SEPULCRE with transfer matrix A=Υ−K ⋅ IN where K is a vector containing the turnover rates for each compartment such that, for compartment j, kj =𝜆 j + ∑N i=1 𝜐 i,j , and IN is the identity matrix. 2.2 | Likelihood calculation As mentioned above, tracer addition data usually comprise time series of the compartment sizes as well as their labelled proportions. We use interchangeably the terms ‘compartments’ and ‘pools’, and the terms ‘tracer’, ‘labelled’ element and ‘marked’ 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. To model compartment sizes and labelled proportions, it is necessary to distinguish between two subpopulations of isotopes: the marked population (usually the heavy isotope) and the unmarked population (usually the light population). Following the notations from López- Sepulcre et al. (2020), we define the state of those subpopulations in the network at time t by two vectors m (t)= ( m(t) 1,…,m(t) N ) and n (t)= ( n(t) 1,…,n(t) N ) giving the marked and unmarked quantities of atoms, respectively, for each compartment. 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 z(t) = m(t) ⊘ x(t) where ⊘ is the elementwise 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 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 x(t) above, the only differences being that y(t) is split into ym(t) and yn(t) (the vectors defining the input of marked and unmarked matter, respectively, for each compartment) and that initial conditions usually differ for m(t) and n(t) . Additionally, for radioactive tracers, the λ rates are adjusted to take into account the radioactive decay rate. For a given set of parameter values, once the system of differential equations is solved and expected trajectories are calculated for x(t) and z(t) , the observed time series for compartment sizes and labelled proportions are modelled as sampling and measurement errors around the expected trajectories. Several error distributions are implemented in isotracer. For example, for the observed size of compartment i at 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 parameter 𝜙 estimated by the model: Other distributions can be a reasonable choice for the generation of observations from the expected trajectories. For example, Normal distributions truncated at 0 or Gamma distributions can be a good approximation if their standard deviations are small and if the labelled proportions are well below 1, and are offered as an option in isotracer. 2.3 | Modelling experimental design variations The mathematical model described above can be adjusted to reflect specific features of a given tracer addition setup. 2.3.1 | The addition regime A core component of modelling a tracer addition is the addition regime. The design of the addition regime is flexible, and additions can be done either at discrete time points (pulses, e.g. Barthelemy et al., 2018; McRoy & Barsdate, 1970; Rønnestad et al., 2000) or continuously over a given duration (continuous intervals or drips, e.g. Collins et al., 2016; de Goeij et al., 2013; Williams et al., 2016). Additions can consist of only marked material, or a mix of marked and unmarked material. Finally, additions may be followed or not by an addition of fully unmarked material (chase, e.g. Bacher et al., 2016; Carbone et al., 2007; Simard et al., 1997). All those regimes can be specified by an appropriate definition of the elements of vector y(t) i . Some source compartments might be considered to be in a steady state on the timescale of the experiment. This can be the case of dissolved nutrients in a stream which 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 (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 ith row of the A matrix to zero and yi(t) to zero: this will result in a constant value for x(t) i equal to the initial condition x( t0 ) i . 2.3.2 | Overenrichment Some sampled compartments can appear overenriched compared to their source compartments. Overenrichment happens when a receiving compartment appears more labelled than its source compartment. This occurs notably in ecosystem studies. Overenrichment should not be possible according to the model, which assumes that the marked material is instantaneously mixed in the receiving compartment pool and that transfer from one compartment does not preferentially affect marked or unmarked material. However, some sampled compartments might actually represent the average of several subcompartments, for example an active compartment involved in material flows and a refractory one which (4) x (t) obs,i∼Truncated Normallower=0 ( 𝜇=x(t) i,𝜎=cv×x(t) i ) (5) z (t) obs,i∼Beta ( 𝛼=𝜙×z(t) i,𝛽=𝜙× ( 1−z(t) i )). | 1123 Methods in Ecology and Evoluon BRUNEAUX ANd LóPEZ- SEPULCRE behaves as if isolated from the rest of the network on the timescale 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 (the active portion), assimilating dissolved nutrients and preferentially grazed by invertebrate consumers, and a refractory portion of sediment and slowly decaying matter (e.g. wood and leaves), not involved in significant nutrient cycling in the span of the experiment. Invertebrate consumers can appear overenriched compared to the average 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, modelling can be adjusted by adding a 𝜋i parameter representing the active fraction of its pool at t0 , and by only taking into account this active subpool when calculating material flows, but by taking into account the inactive, refractory portion when calculating expected total pool sizes and apparent labelled proportions (López- Sepulcre et al., 2020). 3 | PACKAGE IMPLEMENTATION 3.1 | Overview The model exposed above is implemented in isotracer using a Bayesian approach. The Stan program (Carpenter et al., 2017) and its R interface (the rstan package, Stan Development Team, 2020) are used for posterior sampling. Stan uses a Hamiltonian Monte Carlo (HMC) algorithm to sample parameter posteriors with MCMC sampling. The HMC sampling allows for efficient and robust sampling of the posterior, and Stan output provides diagnostics to check that the sampled Markov chains behave appropriately. The Stan model calculates the expected compartment trajectories for a set of parameter values using matrix exponentials to solve the system of differential equations, and the likelihood is calculated by comparing observations to those expected trajectories. Matrix exponentials can be used to solve the system of differential equations in Equation 3 because all the addition regimes currently implemented by isotracer are equivalent to partitioning the experiment timeline into segments over which y(t)=0 . Note that, while the matrix exponential approach is numerically accurate and is the default solver used by isotracer, runtime can become prohibitively long for large models (i.e. with a large number of compartments and/or a large number of unique observation times). A fallback solver is available to the user: it uses a forward Euler scheme with constant step size, but offers no extensive check of numerical accuracy during the numerical solving. This approach allows for greater speed of model evaluation and is expected to perform robustly when an appropriate dt time step and reasonable parameter priors are chosen, but numerical accuracy is not guaranteed. As such, it should be used only for prototyping purposes when exploring the modelling of large datasets, and the matrix exponential solver should always be used for final runs. 3.2 | Interface overview The isotracer package separates the process of modelling tracer additions into two steps: (a) model definition, and (b) model fitting. This is due to the amount of information required to define a network model, namely: network topology, initial conditions, observations and priors. Figure 2 shows an overview of the typical workflow when using isotracer and Table 1 presents the main functions that are used in the model building step. The package functions have been designed for convenient use with the pipe operator %>% provided by the magrittr package and familiar to tidyverse users, but can also be used without the pipe. A new, empty network model is initialized with the new_network- Model() function. This function returns a tibble data frame with zero rows and the supplementary class networkModel, ready to be populated as the user defines the model in subsequent steps. The isotracer package stores its userfacing objects using a tidy table approach (Wickham, 2014) as much as possible to make them easy to inspect. However, the user should use the package functions to modify a model to ensure it results in a valid model. Each row of a networkModel object corresponds to one network replicate. It is often necessary to be able to handle several replicates of a similar network system (e.g. modelling replicate sampling locations in a lake, modelling several transect locations along a stream, sampling several plants in a study of nutrient allocation or modelling the removal of an injected drug in the bloodstream of several patients). All rows (i.e. all network replicates) share the same network topology, and each row can be considered either as a simple replication unit with its own initial conditions and observations, but sharing the parameter values that govern flows with all other replicates, or as a level of a treatment where categorical covariates are taken into account when estimating the model parameters, which can differ across treatments in this case. In other words, defining parameters as same or different among replicates defines whether there are different treatment or factor levels. 3.3 | Model definition 3.3.1 | A minimal model The usage of the functions that define a network model (Table 1) is presented in more details in the case studies and the corresponding appendix below. Here we give an overview of the minimum steps required to define a minimal executable model, and later of the optional steps which allow finetuning model definition. At least five functions have to be used to define a network model (the functions with an asterisk in the last column of Table 1): • new_networkModel() to initialize a new, empty network- Model object. • set_topo() to define the network topology of a network- Model (i.e. the network compartments and their connections). 1124 | Methods in Ecology and Evoluon BRUNEAUX ANd LóPEZ- SEPULCRE • set_init() to provide the initial conditions (total size and labelled proportion) for each compartment. It takes a data frame with columns for compartment identity, size and labelled proportion and accepts supplementary grouping variables used to define replicates. • set_obs() to provide the observations (sampling time, total size, labelled proportion) for each compartment. It accepts a data frame in the same format as set_init(), with an extra column for time. • set_prior() to set prior distributions for individual model parameters. Implemented priors are Uniform, Normal, half- Cauchy, Exponential, Gamma and scaled Beta distribution. All priors are truncated to 0. Additionally, a ‘constant’ prior is available to fix the value of some parameters during MCMC. FIGURE 2 Workflow overview for network modelling with isotracer | 1125 Methods in Ecology and Evoluon BRUNEAUX ANd LóPEZ- SEPULCRE Importantly, while many packages implementing Bayesian methods provide default priors, isotracer does not: the user has to explicitly set priors for all parameters. This is a conscious design choice made to encourage users to make reasoned, explicit choices about their model priors, rather than trusting some default priors that might be woefully inappropriate for the system being modelled. The tutorials presented in appendix and the package documentation provide some guidance for choosing slightly informative priors. 3.3.2 | Additional network properties Depending on the network system that is being modelled, additional network properties can be set using the following functions: • set_steady() to define which compartments should be considered at a steady state. The total size and labelled fraction of compartments at steady state do not change during the calculation of expected compartment trajectories (except if a pulse event is defined for them). This allows modelling, for example, compartments which are constantly replenished (such as dissolved nutrients in flowing stream water) or whose flows are insignificant relative to the timescale of the experiment. • set_split() to define which compartments should be modelled as comprised of an active portion and a refractory portion. • set_half_life() to define the halflife of radioactive isotopic tracers. This allows to take into account the radioactive decay of the tracer in the estimate of flows. By default, isotracer assumes that a stable tracer is used and no halflife is applied in the model. • add_pulse_event() to add pulse events to specific network compartments. A pulse event is defined by an event time and by the quantities of labelled and unlabelled material added to the compartment at that time. In addition to 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 steadystate compartment will result in a new steady state for this compartment, and ‘dripon’ and ‘dripoff’ phases can thus TABLE 1 isotracer functions used to define a network model (valid for version 1.0) Function Role Example Notes Definition Of Network Properties new_networkModel() Create a new, empty network model m <- new_networkModel() * set_topo() Define a network topology m <- set_topo(m, "NH4 - > algae - > grazer") * set_init() Define initial conditions m <- set _init(m, inits, comp = "compartment", size = "biomass", prop = "fraction", group_by = "stream_id") *, a set_obs() Define observations m <- set _obs(m, obs, time = "time_days") *, a set_steady() Set steadystate compartment m <- set _steady(m, "NH4") set_split() Set split compartment m <- set_split(m, "algae") set_half_life() Set a halflife (for radioactive tracer) m <- set_half_life(m, hl = 34.152) add_pulse_event() Set pulse and drip events (see Trinidadian stream case study for examples) b Definition of statistical properties set_prior() Set prior distributions for parameters m <- set_prior(m, normal_p(O, 5), "lambda algae") *, c add_covariates() Define categorical covariates m <- add_covariates(m, upsilon_ NH4_to_algae ~ stream_id) set_size_family() Set distribution for observed pool sizes m <- set_size_family(m, "normal_sd") d set_prop_family() Set distribution for observed proportions m <- set_prop_family(m, "beta_phi") d *Required to build a runnable model. aThe example assumes that both inits and obs have columns named compartment, biomass, fraction and group _ id, and that obs also has a column named time _ days. b(Steady) drips can be defined by applying pulse events to steadystate compartments. cNote that isotracer does not provide any default prior: the user must explicitly set priors for all model parameters. Tutorials in the appendix and package documentation provide guidance for reasonable prior choices. dDefault families are normal distributions for sizes and gamma distributions for proportions. 1126 | Methods in Ecology and Evoluon BRUNEAUX ANd LóPEZ- SEPULCRE be defined using a succession of pulse events to adjust the steady state of a compartment. 3.3.3 | Additional statistical model properties All the functions in the previous section are used to define network properties. The statistical properties of a network model can be defined with the following functions, in addition to the set_prior() function required to define a minimal runnable model: • add_covariates() to define categorical covariates with fixed effects on the estimated parameters. • set_size_family() to define the distribution used to model measurement and sampling error when calculating the likelihood of observed compartment sizes compared to expected trajectories. Implemented distributions are: • set_prop_family() to define the distribution used to model measurement and sampling error when calculating the likelihood of observed labelled fractions compared to expected trajectories. Implemented distributions are: 3.4 | MCMC and posterior analysis Once a network model is defined as explained above, MCMC sampling is performed 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 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 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, the user is expected to use the default output from run_mcmc() (the mcmc.list object). Table 2 provides an overview of the main functions available for postrun processing of mcmc.list objects. The isotracer package provides additional methods for the mcmc.list class to facilitate the manipulation of such objects. For example, a plot method provides a compact way of visualizing MCMC traces and their histograms, and MCMC traces for derived parameters can be easily calculated from the primary parameters returned by run_mcmc() via methods implementing the common mathematical operators ( + , − , × , /, exp, log) for mcmc.list objects. To assess model validity, isotracer provides functions to perform posterior predictive checks. A predict() method is provided for networkModel objects, which takes a model fit (the mcmc. list output from run_mcmc()) and calculates predicted trajectories based on the parameter posteriors. A tidy_dpp() function Normal ( mean, coefficient of variation) Normal (mean, standard deviation) Gamma ( mean, coefficient of variation) Normal (mean, coefficient of variation ) Normal (mean, standard deviation) Beta (mean, precision) TABLE 2 isotracer postrun functions typically used on the mcmc.list output of run _ mcmc() (valid for version 1.0). The examples below assume that MCMC sampling was run on a network model m with for example: fit <- run _ mcmc(m) Function Role Example Notes Visualization and manipulation of MCMC samples plot() Draw traceplots and histograms plot(fit) summary() Get summary of parameter estimates summary(fit) as.matrix Convert MCMC samples to a matrix as.matrix(fit) tidy_mcmc() Convert MCMC samples to a tidy table tidy_mcmc(fit) [, ... ] Extract parameters loss <- fit[, "lambda_algae"] +,- ,*,/,exp(),log() Calculate derived parameters turnover <- 1 /loss Posterior predictive checks predict() Add posterior predictions to a model (a plot() method is implemented) p <- predict(m, fit, probs = 0.9) plot(p) tidy_dpp() Provide tidy data and posterior predictions (useful for use with e.g. the bayesplot package) dpp <- tidy_dpp(m, fit) Network properties (calculated for each MCMC iteration, as a tidy table) tidy_trajectories() Calculate compartment trajectories tidy_trajectories(m, fit) tidy_steady_states() Calculate steadystate compartment sizes tidy_steady_states(m, fit) a tidy_flows() Calculate flows between compartments tidy_flows(m, fit) b aSteady states can only be calculated for compatible networks (e.g. with at least one steadystate source). bFlows can be averaged over the predicted trajectories or calculated for steady states (if the network admits steady states). | 1133 Methods in Ecology and Evoluon BRUNEAUX ANd LóPEZ- SEPULCRE Opinion in Biotechnology, 64, 92– 100. https://doi.org/10.1016/j. copbio.2019.11.003 Bacher, A., Chen, F., & Eisenreich, W.. (2016). Decoding biosynthetic pathways in plants by pulsechase strategies using 13CO2 as a universal tracer. Metabolites, 6(3), 21. https://doi.org/10.3390/metab o6030021 Barneche, D. R., Hulatt, C. J., Dossena, M., Padfield, D., Woodward, G., Trimmer, M., & Yvon- Durocher, G. (2021). Warming impairs trophic transfer efficiency in a longterm field experiment. Nature, 592(7852), 76– 79. h t t p s : / / d o i . o r g / 1 0 . 1 0 3 8 / s 4 1 5 8 6 - 0 2 1 - 0 3 3 5 2 - 2 Barthelemy, H., Stark, S., Michelsen, A., & Olofsson, J. (2018). Urine is an important nitrogen source for plants irrespective of vegetation composition in an Arctic tundra: Insights from a 15N- enriched urea tracer experiment. Journal of Ecology, 106(1), 367– 378. https://doi. org/10.1111/1365- 2745.12820 Cabana, G., & Rasmussen, J. B. (1994). Modelling food chain structure and contaminant bioaccumulation using stable nitrogen isotopes. Nature, 372(6503), 255– 257. https://doi.org/10.1038/372255a0 Carbone, M. S., Czimczik, C. I., McDuffee, K. E., & Trumbore, S. E. (2007). Allocation and residence time of photosynthetic products in a boreal forest using a lowlevel 14C pulsechase labeling technique. Global Change Biology, 13(2), 466– 477. https://doi. org/10.1111/j.1365- 2486.2006.01300.x Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M. A., Guo, J., Li, P., & Riddell, A.. (2017). Stan: A probabilistic programming language. Journal of Statistical Software, 76(1), 1– 32. https://doi.org/10.18637/ jss.v076.i01 Collins, S. M., Thomas, S. A., Heatherly, T., MacNeill, K. L., Leduc, A. O. H. C., López- Sepulcre, A., Lamphere, B. A., El- Sabaawi, R. W., Reznick, D. N., Pringle, C. M., & Flecker, A. S.. (2016). Fish introductions and light modulate food web fluxes in tropical streams: A wholeecosystem experimental approach. Ecology, 97, 3154– 3166. https://doi.org/10.1002/ecy.1530 Crossley, D. A., & Reichle, D. E. (1969). Analysis of transient behavior of radioisotopes in insect food chains. Bioscience, 19(4), 341– 343. https://doi.org/10.2307/1294518 de Goeij, J. M., van Oevelen, D., Vermeij, M. J. A., Osinga, R., Middelburg, J. J., de Goeij, A. F. P. M., & Admiraal, W. (2013). Surviving in a marine desert: The sponge loop retains resources within coral reefs. Science, 342(6154), 108– 110. https://doi.org/10.1126/scien ce.1241981 Dodds, W. K., Collins, S. M., Hamilton, S. K., Tank, J. L., Johnson, S., Webster, J. R., Simon, K. S., Whiles, M. R., Rantala, H. M., McDowell, W. H., Peterson, S. D., Riis, T., Crenshaw, C. L., Thomas, S. A., Kristensen, P. B., Cheever, B. M., Flecker, A. S., Griffiths, N. A., Crowl, T., … Martí, E. (2014). You are not always what we think you eat: Selective assimilation across multiple wholestream isotopic tracer studies. Ecology, 95(10), 2757– 2767. https://doi. org/10.1890/13- 2276.1 Dodds, W. K., Evans- White, M. A., Gerlanc, N. M., Gray, L., Gudder, D. A., Kemp, M. J., López, A. L., Stagliano, D., Strauss, E. A., Tank, J. L., Whiles, M. R., & Wollheim, W. M. (2000). Quantification of the nitrogen cycle in a prairie stream. Ecosystems, 3(6), 574– 589. https:// doi.org/10.1007/s1002 10000050 Freeman, C. J., Thacker, R. W., Baker, D. M., & Fogel, M. L. (2013). Quality or quantity: Is nutrient transfer driven more by symbiont identity and productivity than by symbiont abundance? The ISME Journal, 7(6), 1116– 1125. https://doi.org/10.1038/ismej.2013.7 Fry, B. (2006). Stable isotope ecology. Springer- Verlag. ISBN 9 7 8 – 0 – 3 8 7 - 3 0 5 1 3 - 4 . Gelman, A., Bois, F., & Jiang, J. (1996). Physiological pharmacokinetic analysis using population modeling and informative prior distributions. Journal of the American Statistical Association, 91(436), 1400– 1412. https://doi.org/10.1080/01621 459.1996.10476708 Gillespie, W. R., Zhang, Y., Margossian, C., & Torsten Development Team. (2021). Torsten, a pharmacokinetics/pharmacodynamics library for Stan (in development). https://github.com/metru mrese archg roup/ torsten Goodale, C. L., Fredriksen, G., Weiss, M. S., McCalley, C. K., Sparks, J. P., & Thomas, S. A. (2015). Soil processes drive seasonal variation in retention of 15N tracers in a deciduous forest catchment. Ecology, 96(10), 2653– 2668. Hershey, A. D., & Chase, M. (1952). Independent functions of viral protein and nucleic acid in growth of bacteriophage. The Journal of General Physiology, 36(1), 39– 56. https://doi.org/10.1085/jgp.36.1.39 Kim, I.- Y., Suh, S.- H., Lee, I.- K., & Wolfe, R. R. (2016). Applications of stable, nonradioactive isotope tracers in in vivo human metabolic research. Experimental & Molecular Medicine, 48(1), e203. https://doi. org/10.1038/emm.2015.97 Kulmatiski, A., Beard, K. H., Verweij, R. J. T., & February, E. C. (2010). A depthcontrolled tracer technique measures vertical, horizontal and temporal patterns of water use by trees and grasses in a subtropical savanna. New Phytologist, 188(1), 199– 209. https://doi. org/10.1111/j.1469- 8137.2010.03338.x Li, L., Nelson, C. J., Troesch, J., Castleden, I., Huang, S., & Millar, A. H.. (2018). Data from: Protein degradation rate in Arabidopsis thaliana leaf growth and development. https://doi.org/10.5061/dryad. q3h85 Li, L., Nelson, C. J., Trösch, J., Castleden, I., Huang, S., & Millar, A. H. (2017). Protein degradation rate in Arabidopsis thaliana leaf growth and development. The Plant Cell, 29(2), 207– 228. https://doi. org/10.1105/tpc.16.00768 López- Sepulcre, A., Bruneaux, M., Collins, S. M., el- Sabaawi, R., Flecker, A. S., & Thomas, S. A. (2020). A new method to reconstruct quantitative food webs and nutrient flows from isotope tracer addition experiments. The American Naturalist, 195(6), 964– 985. https://doi. org/10.1086/708546 Lunn, D. J., Best, N., Thomas, A., Wakefield, J., & Spiegelhalter, D. (2002). Bayesian analysis of population PK/PD models: General concepts and software. Journal of Pharmacokinetics and Pharmacodynamics, 29(3), 271– 307. https://doi.org/10.1023/A:10202 06907668 McRoy, C. P., & Barsdate, R. J. (1970). Phosphate absorption in eelgrass. Limnology and Oceanography, 15(1), 6– 13. https://doi.org/10.4319/ lo.1970.15.1.0006 Mill, A. C., Pinnegar, J. K., & Polunin, N. V. C. (2007). Explaining isotope trophicstep fractionation: Why herbivorous fish are different. Functional Ecology, 21(6), 1137– 1145. https://doi. org/10.1111/j.1365- 2435.2007.01330.x Mulholland, P. J., Tank, J. L., Sanzone, D. M., Wollheim, W. M., Peterson, B. J., Webster, J. R., & Meyer, J. L. (2000). Nitrogen cycling in a forest stream determined by a 15N tracer addition. Ecological Monographs, 70(3), 471– 493. https://doi.org/10.2307/2657212 Norman, B. C., Whiles, M. R., Collins, S. M., Flecker, A. S., Hamilton, S. K., Johnson, S. L., Rosi, E. J., Ashkenas, L. R., Bowden, W. B., Crenshaw, C. L., Crowl, T., Dodds, W. K., Hall, R. O., el- Sabaawi, R., Griffiths, N. A., Marti, E., McDowell, W. H., Peterson, S. D., Rantala, H. M., … Webster, J. R. (2017). Drivers of nitrogen transfer in stream food webs across continents. Ecology, 98(12), 3044– 3055. https://doi. org/10.1002/ecy.2009 Riis, T., Dodds, W. K., Kristensen, P. B., & Baisner, A. J. (2012). Nitrogen cycling and dynamics in a macrophyterich stream as determined by a 15N– NH4 + release. Freshwater Biology, 57(8), 1579– 1591. https:// doi.org/10.1111/j.1365- 2427.2012.02819.x Robinson, D. (2001). δ15N as an integrator of the nitrogen cycle. Trends in Ecology & Evolution, 16(3), 153– 162. https://doi.org/10.1016/S0169 - 5 3 4 7 ( 0 0 ) 0 2 0 9 8 - X Rønnestad, I., Dominguez, R. P., & Tanaka, M. (2000). Ontogeny of digestive tract functionality in Japanese flounder, Paralichthys olivaceus studied by in vivo microinjection: pH and assimilation of free amino acids. Fish Physiology and Biochemistry, 220(3), 225– 235. https:// doi.org/10.1023/A:10078 01510056 1134 | Methods in Ecology and Evoluon BRUNEAUX ANd LóPEZ- SEPULCRE Sánchez- Carrillo, S., & Álvarez Cobelas, M. (2017). Stable isotopes as tracers in aquatic ecosystems. Environmental Reviews, 26(1), 69– 81. https://doi.org/10.1139/er- 2017- 0040 Shik, J. Z., Rytter, W., Arnan, X., & Michelsen, A. (2018). Disentangling nutritional pathways linking leafcutter ants and their coevolved fungal symbionts using stable isotopes. Ecology, 99(9), 1999– 2009. https://doi.org/10.1002/ecy.2431 Simard, S. M., Durall, D., & Jones, M. (1997). Carbon allocation and carbon transfer between Betula papyrifera and Pseudotsuga menziesii seedlings using a 13C pulselabeling method. Plant and Soil, 191, 41– 55. https://doi.org/10.1023/A:10042 05727882 Stan Development Team (2020). RStan: The R interface to Stan. R package version 2.21.2. https://cran.rproje ct.org/packa ge=rstan Tank, J. L., Martí, E., Riis, T., Schiller, D., Reisinger, A. J., Dodds, W. K., Whiles, M. R., Ashkenas, L. R., Bowden, W. B., Collins, S. M., Crenshaw, C. L., Crowl, T. A., Griffiths, N. A., Grimm, N. B., Hamilton, S. K., Johnson, S. L., McDowell, W. H., Norman, B. M., Rosi, E. J., … Webster, J. R. (2018). Partitioning assimilatory nitrogen uptake in streams: An analysis of stable isotope tracer additions across continents. Ecological Monographs, 88( 1 ) , 1 2 0 – 1 3 8 . https://doi.org/10.1002/ecm.1280 Weindl, D., Wegner, A., & Hiller, K. (2015). Metabolomewide analysis of stable isotope labeling— Is it worth the effort? Frontiers in Physiology, 6, 1– 3. https://doi.org/10.3389/fphys.2015.00344 Whiles, M. R., Hall, R. O., Dodds, W. K., Verburg, P., Huryn, A. D., Pringle, C. M., Lips, K. R., Kilham, S. S., Colón- Gaud, C., Rugenski, A. T., Peterson, S., & Connelly, S.. (2013). Diseasedriven amphibian declines alter ecosystem processes in a tropical stream. Ecosystems, 16(1), 146– 157. Wickham, H. (2014). Tidy data. Journal of Statistical Software, 59(10), 1– 23. https://doi.org/10.18637/ jss.v059.i10 Williams, C. M., McCue, M. D., Sunny, N. E., Szejner- Sigal, A., Morgan, T. J., Allison, D. B., & Hahn, D. A. (2016). Cold adaptation increases rates of nutrient flow and metabolic plasticity during cold exposure in Drosophila melanogaster. Proceedings of the Royal Society B, 283(1838), 20161317. https://doi.org/10.1098/rspb.2016.1317 Wollheim, W. M., Peterson, B. J., Deegan, L. A., Bahr, M., Hobbie, J. E., Jones, D., Bowden, W. B., Hershey, A. E., Kling, G. W., & Miller, M. C. (1999). A coupled field and modeling approach for the analysis of nitrogen cycling in streams. Journal of the North American Benthological Society, 18(2), 199– 221. https://doi.org/10.2307/1468461 Yuan, J., Fowler, W. U., Kimball, E., Lu, W., & Rabinowitz, J. D. (2006). Kinetic flux profiling of nitrogen assimilation in Escherichia coli. Nature Chemical Biology, 2(10), 529– 530. https://doi.org/10.1038/ nchem bio816 SUPPORTING INFORMATION Additional supporting information may be found in the online version of the article at the publisher’s website. How to cite this article: Bruneaux, M. & López- Sepulcre, A. (2022). isotracer: An R package for the analysis of tracer addition experiments. Methods in Ecology and Evolution, 13, 1119–1134. https://doi.org/10.1111/2041-210X.13822