Example Data

For this example we use one of the samples in the Morris2023 dataset that comes with PoolDilutionR. It consists of data collected at five time points during a methane pool dilution experiment (see ?Morris2023). These values have been converted from ppm to ml, although any single unit would suffice (ie., not a concentration).

Key to columns:

• id - sample identifier
• time_days - time (in days) between samples
• cal12CH4ml - mL’s of ¹²C-methane, calculated from incubation volume and sample ppm
  
• cal13CH4ml - mL’s of ¹³C-methane, calculated from incubation volume and sample ppm
  
• AP_obs - percent of methane that is ¹³C-methane, $\frac{^{13}C}{(^{12}C+^{13}C)} * 100$

library(PoolDilutionR)

OneSample_dat <- subset(Morris2023, subset = id == "71")
print(OneSample_dat)
#>    id  time_days cal12CH4ml cal13CH4ml   AP_obs
#> 26 71 0.00000000  0.6619600 0.00910000 1.356064
#> 27 71 0.02708333  0.3759093 0.00478686 1.257396
#> 28 71 0.06041667  0.2576457 0.00309738 1.187905
#> 29 71 0.08888889  0.2218850 0.00253422 1.129235
#> 30 71 0.11666667  0.2131561 0.00225264 1.045752

The question we are asking this data: What combination of gross fluxes (methane production and consumption) best fits the observed net fluxes implied by these concentration data?

Solving for gross rates

To solve for gross methane production and consumption for this sample, we call pdr_optimize(), which in turn uses R’s optim() to minimize the error in the observed vs predicted pool size and isotopic composition. Errors are weighted by a combination of data quality and explanatory power (ie., cost_fuction()) and predictions are made using pdr_predict(). See Under the hood. for more information.

results <- pdr_optimize(
  time = OneSample_dat$time_days, # time as a vector
  m = OneSample_dat$cal12CH4ml + OneSample_dat$cal13CH4ml, # total pool size
  n = OneSample_dat$cal13CH4ml, # pool size of heavy isotopologue
  P = 0.1, # inital production rate for optim()
  pool = "CH4", # indicates use of default fractionation constants for methane
  m_prec = 0.001, # instrument precision for total pool size, as standard deviation
  ap_prec = 0.01, # instrument precision for atom percent, as standard deviation
)
#> No frac_P provided; looking up from pdr_fractionation table
#> No frac_k provided; looking up from pdr_fractionation table
#> Estimated k0 = 11.873 from n_slope = -11.635

Notice that pdr_optimize() calculated or looked up three parameters that we didn’t provide: frac_P, the fractionation constant for production; frac_k, the fractionation constant for consumption; and k0, the initial value for the first order rate constant of consumption. The first two are taken from default values for methane as provided in the package’s pdr_fractionation dataset, while the latter is estimated from the data, using the slope of the pool size of heavy isotope over time (n_slope).

print(results)
#> $par
#>         P         k 
#>  7.621946 36.495259 
#> 
#> $value
#> [1] 16.62139
#> 
#> $counts
#> function gradient 
#>      102      102 
#> 
#> $convergence
#> [1] 52
#> 
#> $message
#> [1] "ERROR: ABNORMAL_TERMINATION_IN_LNSRCH"
#> 
#> $initial_par
#>        P        k   frac_P   frac_k 
#>  0.10000 11.87255  0.01000  0.98000 
#> 
#> $initial_other
#> $initial_other$method
#> [1] "L-BFGS-B"
#> 
#> $initial_other$lower
#>     P     k 
#> 0e+00 1e-04 
#> 
#> $initial_other$upper
#>   P   k 
#> Inf Inf 
#> 
#> $initial_other$control
#> list()

pdr_optimize() returns a list with the following entries:

• par - the final optimized values of (in this case) production (P) and the first order rate constant of consumption (k)*, as calculated by optim()
• value, counts, convergence, message - these are all returned by optim()
• initial_par - the initial parameters used for the model
• initial_other - other settings passed to optim(); these typically include parameter bounds and optimization method

*From k, the gross consumption rate can be calculated for any pool size of interest.

This list has a great deal of detail, but we can also get a summarized form in a convenient data frame output format:

pdr_optimize_df(
  time = OneSample_dat$time_days,
  m = OneSample_dat$cal12CH4ml + OneSample_dat$cal13CH4ml,
  n = OneSample_dat$cal13CH4ml,
  P = 0.1,
  pool = "CH4",
  m_prec = 0.001,
  ap_prec = 0.01
)
#> No frac_P provided; looking up from pdr_fractionation table
#> No frac_k provided; looking up from pdr_fractionation table
#> Estimated k0 = 11.873 from n_slope = -11.635
#>     par     value initial.P initial.k initial.frac_P initial.frac_k convergence
#> 1     P  7.621946       0.1  11.87255           0.01           0.98          52
#> 1.1   k 36.495259       0.1  11.87255           0.01           0.98          52

Visual Summary of Multisample Output

Now we show a graphical summary of the results from running these samples through pdr_optimize_df() with the same arguments as shown for a single sample above.

samples <- split(Morris2023, Morris2023$id)

all_results <- lapply(samples, FUN = function(samples) {
  pdr_optimize_df(
    time = samples$time_days,
    m = samples$cal12CH4ml + samples$cal13CH4ml,
    n = samples$cal13CH4ml,
    P = 0.1,
    pool = "CH4",
    m_prec = 0.001,
    ap_prec = 0.01,
    quiet = TRUE,
    include_progress = FALSE
  )
})
all_results <- do.call(rbind, all_results)

Below you can see the relative goodness of fit for final optimized P and k values (colored line) in terms of predicted isotopic signature (Atom% 13C) and pool size (Total Methane) for six samples. Grey lines represent values tested on the way to the final optimized fit. Code for reproducing these graphs is not shown, but is available in the publication associated with this package (citation and DOI pending acceptance).

#> Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
#> ℹ Please use `linewidth` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
#> generated.

Pool size predictions (below) are consistently tight to observations, whereas isotopic signature fits (above) have more variation.

Poor fits can result from a variety of issues, including instrument error, low precision, and incorrect assumptions regarding fractionation. To help address the latter, pdr_optimize() allows users to choose up to four parameters (P, k, fractionation of P frac_p, fractionation of k frac_k) that pdr_optimize() will vary for optimizing model fit. Caution should be used in parameter selection, based on prior knowledge of the biological system.

Optimization of Fractionation

Here we tell pdr_optimize() that the production and consumption values are known (and therefore fixed) and ask it to optimize only the fractionation rates:

# In this example, we assume that P and k are known, but that fractionation rates are unknown.

Sample64_dat <- subset(Morris2023, subset = id == "64")

new_fit <- pdr_optimize_df(
  time = Sample64_dat$time_days,
  m = Sample64_dat$cal12CH4ml + Sample64_dat$cal13CH4ml,
  n = Sample64_dat$cal13CH4ml,
  P = 8,
  k = 1,
  params_to_optimize = c("frac_P", "frac_k"),
  m_prec = 0.001,
  ap_prec = 0.01
)
#> No frac_P provided; looking up from pdr_fractionation table
#> No frac_k provided; looking up from pdr_fractionation table

newFitFracP <- new_fit[new_fit$par == "frac_P", ]$value
newFitFrack <- new_fit[new_fit$par == "frac_k", ]$value

# dataframe with new fit (optimizing P_frac and k_frac)
y_new <- pdr_predict(
  time = Sample64_dat$time_days,
  m0 = Sample64_dat$cal12CH4ml[1] + Sample64_dat$cal13CH4ml[1],
  n0 = Sample64_dat$cal13CH4ml[1],
  P = 8,
  k = 1,
  frac_P = newFitFracP,
  frac_k = newFitFrack
)


dat_new <- cbind(Sample64_dat, y_new)
dat_new$oldAPpred <- incdat_out[incdat_out$id == "64", ]$AP_pred


# graph new fit
ggplot(data = dat_new) +
  geom_point(aes(time_days, AP_pred, shape = "New Prediction", color = "New Prediction"),
    size = 3, fill = "#619CFF"
  ) +
  geom_point(aes(time_days, AP_obs, shape = "Observations", color = "Observations"),
    size = 3
  ) +
  geom_point(aes(time_days, oldAPpred, shape = "Old Prediction", color = "Old Prediction"),
    size = 3
  ) +
  scale_shape_manual(
    name = "Sample 64",
    breaks = c("New Prediction", "Observations", "Old Prediction"),
    values = c("New Prediction" = 21, "Observations" = 16, "Old Prediction" = 1)
  ) +
  scale_color_manual(
    name = "Sample 64",
    breaks = c("New Prediction", "Observations", "Old Prediction"),
    values = c("New Prediction" = "black", "Observations" = "#619CFF", "Old Prediction" = "#619CFF")
  ) +
  ylab("Atom%13C\n")

print(new_fit)
#>        par      value initial.P initial.k initial.frac_P initial.frac_k
#> 1   frac_P 0.01292351         8         1           0.01           0.98
#> 1.1 frac_k 0.97710815         8         1           0.01           0.98
#>     convergence
#> 1             0
#> 1.1           0

(As before, we are using pdr_optimize_df() for tidy data frame output.)

We can see that

• the values of P and k were specified
• initial values for frac_P and frac_k were looked up from the package’s pdr_fractionation table based on the pool (here, methane)
• as requested, the optimizer only optimized frac_P and frac_k

The help page for pdr_optimize() provides other examples of this.

Prediction

The equations in pdr_predict() come from the following equations in von Fischer and Hedin 2002.

Equation 5 describes the change in pool size over time in terms of P, gross production, and k, the first order rate constant of consumption.

$$m_t = \frac{P}{k} - (\frac{P}{k} - m_0) * e^{(-kt)}$$

Where:
m_t is the total pool size at time t.
m₀ is the total pool size at time zero.


Equation 9 tracks the change in the heavy isotopologue over time in terms of consumption and its fractionation constant.

n_t = n₀ * e^(−kαt)

Where:

n_t is the pool size of the heavy isotopologue at time t.
n₀ is the pool size of the heavy isotopologue at time zero.

And:
$$ \alpha = \frac{k_{(\,^{13}C)\,}}{k_{(\,^{12}C)\,}} $$
Where:
k_( ¹³C) is the first-order rate constant for consumption of ¹³CH₄
k_( ¹²C) is the first-order rate constant for consumption of ¹²CH₄.


In PoolDilutionR, we expand this equation to allow for the production of heavy molecules from authocthonous sources.


Equation 9, adjusted
$$n_t = \frac{p_{frac}}{k_{frac}} - ( \frac{p_{frac}}{k_{frac}} - 0) * e^{-k_{frac}t}$$

Where k_frac is equivalent to k * α above, for whatever heavy isotope is applicable.


The isotopic composition of the pool over time is described in Equation 10.

$$ AP_t = \frac {n_t}{m_t} + AP_p $$

Which we can now simplify to:

Equation 10, adjusted $$AP_t = (\frac{n_t}{m_t}) * 100$$

Due to production of heavy isotopologue now being accounted for in our adjusted Equation 9.


Combined this looks like:

$$ AP_t = \left( \frac{\frac{p_{frac}}{k_{frac}} - ( \frac{p_{frac}}{k_{frac}} - 0) * e^{-k_{frac}t}} {\frac{P}{k} - (\frac{P}{k} - m_0) * e^{(-kt)}} \right) * 100 $$

Cost Function

The pdr_cost() provides feedback to pdr_predict() on the quality of each iteration of fitted rates and/or fractionation constants. The sum of errors are weighted by the standard deviation of the observations, as well as a scaling factor, (N).


Equation 14:

$$E = \left(\sum_{t=1}^j\frac {AP_{obs}(t) - AP_{pred}(t)}{SD_{obs-AP}}\right) * N_{ap} + \left(\sum_{t=1}^j\frac {m_{obs}(t) - m_{pred}(t)}{SD_{obs-m}}\right) * N_{m}$$

Where:

SD_{obs − AP} is the standard deviation among all observations of atom percent for a single sample
SD_{obs − m} is the standard deviation among all observations of total pool size for a single sample

And:

$$N_x = \frac{SD_{x_{observed}}}{SD_{x_{precision}}}$$
Where:


x is either atom percent (ap) or total pool size (m)
SD_xprecision is the instrument precision for that variable as standard deviation (ie., standard precision)


Users have the ability to replace the default cost function (pdr_cost()) with their own, if desired.

Fractionation

Fractionation describes the tendency of heavier isotopes and isotopologues to resist chemical transformation, whether these be phase changes or chemical reactions, whether spontaneous or enzyme-mediated. There are two major types of fractionation. Equilibrium fractionation occurs when phase changes favor the heavier isotopologue staying a lower energy state (Druhan, Winnick, and Thullner 2019, Urey 1947). A typical example would be the relative enrichment of the ocean in heavy water H₂¹⁸O due to preferential evaporation of light water H₂¹⁶O. The second major type of fractionation is kinetic, and is classically associated enzyme selectivity. This is what drives the distinction in ¹³C signatures between C3 and C4 plants (O’Leary 1981).

Because our knowledge of earth system processes and enzyme diversity is rapidly expanding, additional fractionation constants will be added to pdr_fraction as part of future package versions.