In smwindecker/deconvolve: Mixture Model for the Deconvolution of Thermal Decay Curves

knitr::opts_chunk$set(echo = TRUE)

Background

Plant cell wall biomass is composed of a range of different types of carbon. Ecologists can use these proportions of biomass carbon to make generalisations about species' effects on ecosystem processes such as litter decomposition, which links plant biomass to the global carbon cycle [@cornwell2008]. Biofuels researchers use proportions of biomass carbon types to estimate kinetic decay parameters of species' tissue. Traditional methods for calculation of lignocellulosic biomass involve wet chemistry assays for carbon component analysis, which are time-consuming and adversely impact the environment through use of sulfuric acid and acetic anhydride, among other chemicals. Thermogravimetric analysis (TGA) is an alternative method, already in use among biofuels researchers, to approximate these carbon compounds from mass loss data obtained by heating a biomass sample (in a $N_2$ environment, termed pyrolysis).

Mass loss during pyrolysis can be regarded as the sum of the degradation of the main components of the sample --- hemicelluloses, cellulose, and lignin [@hu2016; @cheng2015; @perejon2011; @orfao2001; @mullerhagedorn2007]. Therefore, the multi-peaked rate of mass loss curve can be mathematically separated into constituent parts with a mixture model --- termed deconvolution. By integrating under these individual peaks, we can estimate the contribution of each component to the overall initial sample. This method has been validated by studies comparing estimates to experimental measurements [@yang2006]. Plant trait measurement is guided by the principle of standardised, widely reproducible methods [@perezharguindeguy2013], but much of the published literature use commercial software to deconvolve the rate of mass loss curve, for example OriginPro [@chen2017], PeakFit [@perejon2011], Fityk [@perejon2011], or Datafit [@cheng2015]. Software accessibility and reproducibility may be one reason why thermogravimetric analysis for carbon type estimation has not yet been widely adopted by functional ecologists, despite being previously used to identify plant species' recalcitrance, for example in marine and coastal macrophytes [@trevathantackett2015], and eucaplyts [@orfao2001]. In addition, proprietary software inhibits researchers already conducting this type of analysis from making their analyses independently reproducible.

The mixchar package provides an open-source set of functions to perform this deconvolution. Although the nonlinear mixture model used for peak separation at the core of this package could be used for many different purposes, mixchar provides specific guidelines related to thermal decay curve analysis and carbon component estimation.

---

Litter collection and preparation

devtools::install_github('ropenscilabs/icon')

r icon::fa('leaf', size = 2) Collect litter

We developed and tested the functions in this package using the thermogravimetric decay data for litter of 29 different plant species. Two species from this set are available as datasets in the package --- the freshwater reed Juncus amabilis (accessed as juncus) and the freshwater fern Marsilea drumondii (accessed as marsilea).

r icon::fa('burn', size = 2) Dry litter

To ensure component estimates are an accurate representation of the original composition of the litter, it is important to dry samples as quickly as possible to prevent decomposition. Plant litter collected for this analysis was placed in moist plastic bags and stored in dark coolers until transported to the lab, and then moved to a dark, refrigerated room. We dried the fresh litter at 60 °C for 72 hours.

r icon::fa('vials', size = 2) Grind litter

Dry litter must be ground in order to be used in thermogravimetric analysis. We ground litter to < $40 \mu m$ using a Retsch Centrifugal Mill ZM200.

---

Thermogravimetric analysis

r icon::fa('box', size = 2) Pyrolyse samples

We pyrolysed 10--20 mg subsamples of dry, ground litter in an N$_2$ environment from 30--800 °C at a temperature ramp of 10 °C/min using a Netzsch TGA-FTIR thermogravimetric analyser (Department of Biomedical Engineering, University of Melbourne). The resulting data is mass against temperature (Fig. 1).

library(mixchar)
tmp <- process(juncus, init_mass = 18.96,
               temp = 'temp_C', mass_loss = 'mass_loss')
plot(tmp, plot_type = 'mass')

---

Deconvolution

1. Load data

The first step is to load the data derived from the thermogravimetric analysis into R. This step will vary depending on the format of your exported data. Two of many options for reading data are read.csv() for .csv files and read.table() for .txt files. Both of these functions also have a skip argument, in which you can tell R to begin reading data from a certain line, useful if the exported data has several lines of metadata at the top. You may also then have to assign column names. For example:

``` {r, eval = FALSE} my_data <- read.csv('your_file_path_here.csv', header = FALSE, skip = 15) colnames(c('temperature_C', 'mass_mg'))

More information about reading data into R can be found [here](http://www.sthda.com/english/wiki/reading-data-from-txt-csv-files-r-base-functions). 

This package was developed using data from a Netzsch TGA-FTIR thermogravimetric analyser. We have included two example datasets in this package, called `juncus` and `marsilea`, that we will use for this demonstration. 

``` {r}
library(mixchar)
head(juncus)

---

2. Calculate rate of mass loss

After we've loaded our data, we need to calculate the rate of mass loss across temperature by taking the derivative. We use the process function to take the derivative, resulting in rate of mass loss over temperature data. The process function needs the dataset, the initial mass of sample, the name of the temperature data column, and the name of your mass column (mg). Mass might be recorded in your exported dataset in different ways:

• If you have mass loss of sample, such as in the juncus dataset, you can specify the column name using the mass_loss argument.
• If you have mass remaining of sample, you specify the column name using the mass argument.
• You only need to specify one of these columns for the process function. But make sure the argument matches the kind of data!

The function defaults to temperature data in Celsius, but you can also modify to indicate the data is provided in Kelvin, by specifying the argument temp_units = 'K'.

deriv_juncus <- process(juncus,                          # dataframe name
                        init_mass = 18.96,
                        temp = 'temp_C',                 # column name for temperature
                        mass_loss = 'mass_loss')         # column name for mass loss data
deriv_juncus

If you plot() the output of the process function, you will get two curves: the mass of sample across time and the rate of mass loss curve. If you're only interested in one plot, you can specify plot_type = 'mass' or plot_type = 'rate'. You can also specify a size factor using cex. The rate of mass loss curve is a multi-peaked curve (Fig. 2) encompassing three main phases [@orfao2001]:

1. A short period with a pronounced peak of moisture evolution, up until approximately 120 °C.
2. A wide mid-range of high mass loss, caused by devolatilisation of primary biomass carbon components, between approximately 120--650 °C.
3. A final period of little mass loss when carbonaceous material associated with the inorganic fraction decomposes, after approximately 650 °C.

library(mixchar)
tmp <- process(juncus, init_mass = 18.96,
               temp = 'temp_C', mass_loss = 'mass_loss')
plot(tmp, plot_type = 'rate')
segments(x0 = 40, y0 = 0.0018, x1 = 40, y1 = 0.002, lwd = 3, col = 'darkgrey')
segments(x0 = 40, y0 = 0.002, x1 = 120, y1 = 0.002, lwd = 3, col = 'darkgrey')
segments(x0 = 120, y0 = 0.0018, x1 = 120, y1 = 0.002, lwd = 3, col = 'darkgrey')
text(x = ((120-40)/2+40), y = 0.0024, '1', cex = 2, col = 'darkgrey')

segments(x0 = 120, y0 = 0.0078, x1 = 120, y1 = 0.008, lwd = 3, col = 'darkgrey')
segments(x0 = 120, y0 = 0.008, x1 = 650, y1 = 0.008, lwd = 3, col = 'darkgrey')
segments(x0 = 650, y0 = 0.0078, x1 = 650, y1 = 0.008, lwd = 3, col = 'darkgrey')
text(x = 400, y = 0.0076, '2', cex = 2, col = 'darkgrey')

segments(x0 = 650, y0 = 0.0008, x1 = 650, y1 = 0.001, lwd = 3, col = 'darkgrey')
segments(x0 = 650, y0 = 0.001, x1 = 790, y1 = 0.001, lwd = 3, col = 'darkgrey')
segments(x0 = 790, y0 = 0.0008, x1 = 790, y1 = 0.001, lwd = 3, col = 'darkgrey')
text(x = ((790-650)/2+650), y = 0.0014, '3', cex = 2, col = 'darkgrey')

In the next step we will crop the data to include only phase 2, so visualising your own data is important to check that the default temperature bounds will be suitable, in case for example the dehydration phase extends past 120 $^{\circ}$C or ends earlier. Comparing the plots for the two species we can see similarities in the shape and location of the peaks of the overall rate of mass loss curve, but also subtle differences. It is these characteristics we will tease apart using the nonlinear mixture model in the next step.

---

3. deconvolve data

Crop mass loss data

The decovolve function takes care of modelling the rate of mass loss data with the nonlinear mixture model. Since the overall DTG curve thus represents the loss of extractives, water, inorganic matter, and volatiles in addition to our components of interest [@hu2016], we isolate mass loss from our primary biomass components by cropping the DTG data to Phase 2 (roughly 120--650 °C). The default temperature bounds are 120 $^{\circ}$C and 700 $^{\circ}$C, but can be modified with the lower_temp and upper_temp arguments.

Non-linear mixture model

Biomass components decompose relatively independently because they do not interact much during thermal volatilisation [@yang2006]. Therefore, the cropped DTG curve can be mathematically deconvolved into constituent parts using a mixture model. The derivative rate of mass loss equation ($\frac{dm}{dT}$) can be expressed as the sum of $n$ independent reactions, as follows [@orfao2001]:

\begin{align} -\frac{dm}{dT} &= \sum\limits_{i=1}^n c_i\frac{d\alpha_{i}}{dT} \label{eqn:mixture_model} \ m &= \frac{M_T}{M_0} \label{eqn:fraction} \ c_i &= M_{0i} - M_{\infty i} \label{eqn:decayed_mass} \ \alpha_i &= \frac{M_{0i} - M_{Ti}}{M_{0i} - M_{\infty i}} \label{eqn:alpha} \end{align}

where mass ($m$) is expressed as a fraction of mass at temperature $T$ ($M_T$) of the initial sample mass ($M_0$), $c_i$ is the mass of component $i$ that is decayed, and the mass loss curve of each individual component ($\frac{d\alpha_{i}}{dT}$) is the derivative of $\alpha_i$, the conversion of mass at a given temperature ($M_{Ti}$), from the initial ($M_{0i}$), given total mass lost between the initial and final ($M_{\infty i}$) temperature for each curve.

Although most of our results can be described with only $n = 3$ peaks, corresponding to a single curve each of hemicelluose, cellulose, and lignin, some species yield a second hemicellulose peak at a lower temperature, resulting in $n = 4$ independent curves. This is because the soluble carbohydrates in plant tissue can take many forms, including xylan, amylose, etc., which apparently degrade at different temperatures [see also @chen2017; @mullerhagedorn2007]. deconvolve will decide whether three or four peaks are best using an internal function that determines if there is a peak below 220 $^{\circ}$C. Upon inspection of your curve you can override this by modifying the n_peaks argument.

In order to fit the mixture model, we must determine the shape of the individual curves ($\frac{d\alpha_{i}}{dT}$) that are summed to produce it. Many different functions have been proposed: the asymmetric bi-Gaussian [@sun2015], logistic [@barbadillo2007], Weibull [@cai2007], asymmetric double sigmoidal [@chen2017], and the Fraser-Suzuki function [@perejon2011; @hu2016]. Comparisons of several techniques [@svoboda2013; @perejon2011; @cheng2015] found that the Fraser-Suzuki function best fit these kinetic curves, since it allows for asymmetry. We therefore use the Fraser-Suzuki function to describe the rate expression of a single curve as follows:

\begin{gather}\label{eqn:fs_function} \frac{d\alpha_i}{dT} = h_i\ exp\bigg{-\frac{ln2}{s_i^2}\Big[ln\Big(1 + 2s_i \frac{T - p_i}{w_i}\Big)\Big]^2\bigg} \end{gather}

where T is temperature (°C), and the parameters $h_i$, $s_i$, $p_i$, and $w_i$ are height, skew, position, and width of curve $i$, respectively. In total, our model estimates 12 or 16 parameters, one for each parameter for either 3 or 4 primary components. You can play around with changing individual parameters for each curve, and check out the effect on the overall decay curve, using the shiny app.

deconvolve() uses non-linear optimisation with residual sum of squares to fit the rate expression [as in @cheng2015]. Starting values were selected based on curves depicted in the literature [@mullerhagedorn2007] and from the results of running an identical deconvolution on pure cellulose and lignin samples. Hemicelluloses decay in a reasonably narrow band beginning at a lower temperature [@mullerhagedorn2007], so we used 270 for position and 50 for width. Linear cellulose crystals decay at a higher temperature, but decay more rapidly after peak temperatures are reached, so starting position was set to 310 and width to 30. Lignin typically decays beginning at a high temperature and over a wide interval [@mullerhagedorn2007], so position and width began at 410 and 200, respectively. These starting values can also be modified in the deconvolve() function with the start_vec, lower_vec, and upper_vec arguments. These initial starting values are optimised before model fitting using the NLOPTR_LN_BOBYQA algorithm [@bobyqa] within the nloptr [@nloptr] package. Finally, optimised starting values are fit to the non-linear mixture model using nlsLM() in the minpack.lm [@minpack.lm] package.

The function also has built in starting values for the nonlinear optimisation. These values were tested on litter from 30 plant species, encompassing herbaceous, graminoid, as well as woody species. However, if they do not work for your sample, or you would like to play with the effect of changing them, you can do so with the start_vec, lower_vec, and upper_vec arguments.

output_juncus <- deconvolve(deriv_juncus)
output_juncus

deriv_marsilea <- process(marsilea, 
                          init_mass = 15.29,
                          temp = 'temp_C', 
                          mass_loss = 'mass_loss')

# here's an example of specifying your starting vector. 
# the order of values is height, skew, position, then width

my_starting_vec <- c(height_0 = 0.002, skew_0 = -0.15, position_0 = 210, width_0 = 50,
                     height_1 = 0.003, skew_1 = -0.15, position_1 = 250, width_1 = 50,
                     height_2 = 0.006, skew_2 = -0.15, position_2 = 320, width_2 = 30,
                     height_3 = 0.001, skew_3 = -0.15, position_3 = 390, width_3 = 200)

output_marsilea <- deconvolve(deriv_marsilea, n_peaks = 4, start_vec = my_starting_vec)

output_marsilea

Component weights

Once overall curve parameters are fit, we can pass each components's parameter estimates to a single Fraser-Suzuki function, and calculate the weight of the component in the overall sample by integrating the area under the peak. To estimate the uncertainty of the weight predictions, deconvolve() will calculate the 95% interval of the weight estimates across a random sample of parameter estimates, drawn in proportion to their likelihood. We assume a truncated multivariate normal distribution, since the parameters are constrained to positive values, using the modelling package tmvtnorm [@tmvtnorm].

\begin{gather}\label{eqn:integration} \alpha_i = \int_{120}^{650} h_i\ exp\bigg{-\frac{ln2}{s_i^2}\Big[ln\Big(1 + 2s_i \frac{T - p_i}{w_i}\Big)\Big]^2\bigg} dT \end{gather}

We interpret that the peak located around 250--270 °C corresponds to primary hemicelluloses, around 310--330 °C to cellulose, and around 330--350 °C to lignin. If present, the fourth peak located below 200 °C corresponds to the most simple hemicelluloses. The second dataset included in the package, marsilea, is an example of a four-peak deconvolution.

5. Examine outputs

deconvolve results in a few different outputs that you can retrieve with accessor functions.

• rate_data() will show you the processed dataset that results from the process function, useful if you want to play around with other modelling approaches or plotting options:

juncus_rate <- rate_data(output_juncus)
head(juncus_rate)

*temp_bounds() will print the temperature values at which the data were cropped for analysis:

temp_bounds(output_juncus)

• model_fit() will show you the output of the mixture model and estimated parameter values for height, skew, position, and width for each peak. Peak 1 is hemicellulose, peak 2 is cellulose, and peak 3 is lignin. If present, the optional fourth peak located at the lowest temperature interval will be listed as peak 0:

model_fit(output_juncus)

model_fit(output_marsilea)

• component_weights() will display the mean, upper, and lower confidence intervals for each estimated component:

component_weights(output_juncus)

component_weights(output_marsilea)

---

6. Plot deconvolved curves

The default plotting function for the output of the deconvolve function shows you your raw mass data, the estimated full curve from the mixture model, and also plots the individual component peaks using their parameter estimates from the model. The default plot is in black and white, but you can plot a colour version that uses colour-blind friendly viridis colours by specifying bw = FALSE.

Take a look at your plot. Do the estimates seem reasonable?

plot(output_juncus)

plot(output_marsilea, bw = FALSE)

If the estimated peaks do not match your data well, it may be that your sample is a bit different than those used to develop the package. That's not a problem, but you may need to try the deconvolution again with different starting values. If you aren't managing to produce reasonable estimates, get in touch with us via the issues page.

If you want to modify the aesthetics of this plot for your own work, then you can access the parameter estimates as follows:

juncus_parameters <- model_parameters(output_juncus)

and use fs_function() to plot individual component peaks, and fs_mixture function to plot the overall mixture curve.

---

References

smwindecker/deconvolve documentation built on Sept. 20, 2023, 1:49 a.m.