Advanced example

This example emphasizes two scenarii not covered in the simple example:

  • using charge passed over time data, which can be acquired when performing electrochemical measurements (chronocoulometry), and that can help improve the accuracy of the time constants obtained from the fit
  • fitting when the concentration evolution over time experimental data only covers some of the species considered in the model

Kinetic model

The kinetics model we use leads to the following set of differential equations:

\[\rm \frac{d[HMF]}{dt} = -(k_{\textbf{1}1} + k_{\textbf{1}2} + k_{D\textbf{1}})[HMF]\]
\[\rm \frac{d[DFF]}{dt} = k_{\textbf{1}1}[HMF] - (k_{\textbf{2}1} + k_{D\textbf{2}1})[DFF]\]
\[\rm \frac{d[HMFCA]}{dt} = k_{\textbf{1}2}[HMF] - (k_{\textbf{2}2} + k_{D\textbf{2}2})[HMFCA]\]
\[\rm \frac{d[FFCA]}{dt} = k_{\textbf{2}1}[DFF] + k_{\textbf{2}2}[HMFCA] - (k_{\textbf{3}} + k_{D\textbf{3}})[FFCA]\]
\[\rm \frac{d[FDCA]}{dt} = k_{\textbf{3}}[FFCA] - k_{D\textbf{4}}[FDCA]\]
\[\begin{split}\rm \frac{d[D]}{dt} = k_{D\textbf{1}}[HMF] + k_{D\textbf{2}1}[DFF] + k_{D\textbf{2}2}[HMFCA]\\ \rm + k_{D\textbf{3}}[FFCA] + k_{D\textbf{4}}[FDCA] - k_{D^*}[D]\end{split}\]
\[\rm \frac{d[D^*]}{dt} = k_{D^*}[D]\]

In order to write a function describing this model we need to define all the species whose concentrations evolution over time are considered in this study. In that case 5 species concentrations were tracked by HPLC measurements (HMF, DFF, HMFCA, FFCA, FDCA). The degradation products (D and D*) concentrations could not be tracked. For the untracked species we need to consider 10 additional concentrations. Indeed, in the differential equations mentionned above only the total concentration evolution of D and D* are considered. However, these total D and D* concentrations need to be splitted here into the individual concentrations coming from each individual degradation product formation reactions. This is required because humins species formed through different reaction pathways will involve a different amount of charge passed in our electrochemical measurement. With these individual concentrations we are able to calculate the amount of charge passed over time allowing us to include the corresponding experimental data in the fit.

We define the species names considered in this experiment as follow:

species_tracked = [
    "HMF", "DFF", "HMFCA", "FFCA", "FDCA"
]
species_untracked = [
    "D_HMF",  "D_DFF",  "D_HMFCA",  "D_FFCA",  "D_FDCA",
    "Dx_HMF", "Dx_DFF", "Dx_HMFCA", "Dx_FFCA", "Dx_FDCA"
]
species = species_tracked.copy()
species.extend(species_untracked)

HMF, DFF, HMFCA, FFCA and FDCA are the species whose concentrations are tracked experimentally. The D and D* coming from these species are denoted with the same names prepended respectively with a “D_” and a “Dx_”. To sum up we can now write:

\[\begin{split}\rm [D] = [D\_HMF] + [D\_DFF] + [D\_HMFCA]\\ \rm + [D\_FFCA] + [D\_FDCA]\end{split}\]
\[\begin{split}\rm [D^*] = [Dx\_HMF] + [Dx\_DFF] + [Dx\_HMFCA]\\ \rm + [Dx\_FFCA] + [Dx\_FDCA]\end{split}\]

Therefore, the derivatives presented above for the concentration of humins are rewritten in the following way:

\[\rm \frac{d[D\_HMF]}{dt} = k_{D\textbf{1}}[HMF] - k_{D^*}[D\_HMF]\]
\[\rm \frac{d[D\_DFF]}{dt} = k_{D\textbf{1}}[DFF] - k_{D^*}[D\_DFF]\]
\[\rm \frac{d[D\_HMFCA]}{dt} = k_{D\textbf{1}}[HMFCA] - k_{D^*}[D\_HMFCA]\]
\[\rm \frac{d[D\_FFCA]}{dt} = k_{D\textbf{1}}[FFCA] - k_{D^*}[D\_FFCA]\]
\[\rm \frac{d[D\_FDCA]}{dt} = k_{D\textbf{1}}[FDCA] - k_{D^*}[D\_FDCA]\]

The derivatives for the concentration of D* are rewritten in this way:

\[\rm \frac{d[Dx\_HMF]}{dt} = k_{D^*}[D\_HMF]\]
\[\rm \frac{d[Dx\_DFF]}{dt} = k_{D^*}[D\_DFF]\]
\[\rm \frac{d[Dx\_HMFCA]}{dt} = k_{D^*}[D\_HMFCA]\]
\[\rm \frac{d[Dx\_FFCA]}{dt} = k_{D^*}[D\_FFCA]\]
\[\rm \frac{d[Dx\_FDCA]}{dt} = k_{D^*}[D\_FDCA]\]

We can now write the derivatives function derived from these differential equations:

def derivatives(y, t, p):

    """Calculates the derivatives of the concentrations at t

    Used scipy.integrate.odeint to numerically solve the differential
    equations in a given time range.

    Lists ("y" and "dy") used by scipy.integrate.odeint are converted
    to dictionaries ("c" and "dc") in order to make the differentials
    easier to write and read for humans.

    Parameters:
        y (list): concentration values at t
        t (float): time value where the derivatives are calculated
        p (dict): dictionary containing the parameters used to
        calculate the derivatives e.g. time constants
    """

    c = {s:y[i] for i, s in enumerate(species)}

    dc = dict()

    dc["HMF"]      =  - (p["k11"] + p["k12"] + p["kD1"])*c["HMF"]
    dc["DFF"]      = p["k11"]*c["HMF"]                       - (p["k21"] + p["kD21"])*c["DFF"]
    dc["HMFCA"]    = p["k12"]*c["HMF"]                       - (p["k22"] + p["kD22"])*c["HMFCA"]
    dc["FFCA"]     = p["k21"]*c["DFF"] + p["k22"]*c["HMFCA"] - (p["k3"] + p["kD3"])*c["FFCA"]
    dc["FDCA"]     = p["k3"]*c["FFCA"]                       - p["kD4"]*c["FDCA"]

    dc["D_HMF"]    = p["kD1"]*c["HMF"]    - p["kDx"]*c["D_HMF"]
    dc["D_DFF"]    = p["kD21"]*c["DFF"]   - p["kDx"]*c["D_DFF"]
    dc["D_HMFCA"]  = p["kD22"]*c["HMFCA"] - p["kDx"]*c["D_HMFCA"]
    dc["D_FFCA"]   = p["kD3"]*c["FFCA"]   - p["kDx"]*c["D_FFCA"]
    dc["D_FDCA"]   = p["kD4"]*c["FDCA"]   - p["kDx"]*c["D_FDCA"]

    dc["Dx_HMF"]   = p["kDx"]*c["D_HMF"]
    dc["Dx_DFF"]   = p["kDx"]*c["D_DFF"]
    dc["Dx_HMFCA"] = p["kDx"]*c["D_HMFCA"]
    dc["Dx_FFCA"]  = p["kDx"]*c["D_FFCA"]
    dc["Dx_FDCA"]  = p["kDx"]*c["D_FDCA"]

    dy = [dc[name] for name in species]

    return dy

We can then convert the concentrations evolution over time into the charge passed over time using this equation:

\[\rm Q = e N_A V \sum_i n_i C_i\]

With e the electron charge in Coulombs, \(\rm N_{A}\) the Avogadro number, V the volume of solution, \(\rm n_{i}\) the number of charge passed to make one molecule of i, and \(\rm C_{i}\) the concentration of species i. This equation is used to define a c_to_q() function that convert a list of all concentrations evolution over time c into a list of charge passed over time q as follow:

import numpy as np
import scipy.constants as constants

def c_to_q(c):

    """Calculates the charge passed from the concentrations vs time

    Parameters:
        c (numpy.ndarray):
            Concentrations evolution over time, axis 0 is time,
            axis 1 is species.
    """

    # the concentrations are monitored in micromoles/L
    # so we convert them to moles/L
    c *= 1e-6

    # calculate the product ni*Ci for each species
    ni_Ci = list()
    for i, s in enumerate(species):
        ni_Ci.append(2*(i%5 + int(i/5))*c[:,i])

    # calculate the sum of ni*Ci for all species
    sum_ni_Ci = np.sum(ni_Ci, axis = 0)

    # solution volume in L
    V = 100e-3

    # charge passed in C
    q = constants.e*constants.N_A*V*sum_ni_Ci

    return q

Now that the model is defined we can load the raw data and fit it.

Load concentrations and charge passed evolution over time

Here three .csv files are used to get the measured evolution of HMF, DFF, HMFCA, FFCA and FDCA concentrations over time. Three additional .csv files are used for the charge passed over time. These files are loaded in an object of class data.Dataset. Yon can consult the recomendations for the .csv files formatting in the data.Dataset.load_c() and the data.Dataset.load_q() methods documentation.

from chemical_kinetics import data

folders = [f"data/run{i}/" for i in range(1,4)]
files_c = [f"{folder}Reaction Monitoring.csv" for folder in folders]
files_q = [f"{folder}Charge Passed.csv" for folder in folders]

ds = data.Dataset(
    files_c = files_c,
    files_q = files_q,
    t_label = "Time [h]",
    c_label = r"Concentration [$\rm\mu$M]",
    q_label = "Charge passed [C]"
)

Fitting

The only parameters to be fitted are the time constants. The initial concentrations are fixed to the initial values recorded by HPLC for the tracked species. For the untracked species the initial concentrations are fixed to 0.

We first define the time constants parameters:

parameter_args = dict(value = 0.05, min = 0)
parameter_names = [
    "k11","k12","k21","k22","k3",
    "kD1","kD21","kD22","kD3","kD4",
    "kDx"
    ]
parameters = {name: parameter_args for name in parameter_names}

Then we define the initial concentrations parameter for the tracked species:

c0 = {name: dict(vary = False) for name in species_tracked}

And finally, we define the initial concentrations for the untracked species (note that we need here to define an ordered dictionary):

from collections import OrderedDict
c0_untracked = OrderedDict({
    name: dict(value = 0, vary = False) for name in species_untracked
})

With the data loaded, the model defined and the parameters / initial concentrations initialized we can proceed with the fit, using the function fit.fit_dataset():

from chemical_kinetics import fit

fit.fit_dataset(
    dataset = ds,
    derivatives = derivatives,
    c_to_q = c_to_q,
    parameters = parameters,
    c0 = c0,
    c0_untracked = c0_untracked
)

Fit succeeded.

Fit results

The fit results are summarized in the table below. Note that all the initial concentrations (parameters starting with the string c0_) are fixed so only the ks values are varying parameters in this fit. The parameter values resulting from the fit can be printed using the function fit.print_result():

fit.print_result(ds)
name value stderr stderr/value % init. val. vary min max
0 k11 0.00554 0.000195 3.51 0.05 True 0 inf
1 k12 0.00182 5.24e-05 2.89 0.05 True 0 inf
2 k21 0.0649 0.00411 6.33 0.05 True 0 inf
3 k22 0.038 0.00853 22.4 0.05 True 0 inf
4 k3 0.0073 0.000393 5.38 0.05 True 0 inf
5 kD1 0.0302 0.000367 1.22 0.05 True 0 inf
6 kD21 0.0435 0.00708 16.3 0.05 True 0 inf
7 kD22 2.93e-05 0.00839 2.86e+04 0.05 True 0 inf
8 kD3 0.0545 0.00431 7.9 0.05 True 0 inf
9 kD4 0.0727 0.00592 8.13 0.05 True 0 inf
10 kDx 0.00173 0.000275 15.9 0.05 True 0 inf
11 c0_HMF 4.86e+03 0 0 4.86e+03 False -inf inf
12 c0_DFF 4.57 0 0 4.57 False -inf inf
13 c0_HMFCA 0.453 0 0 0.453 False -inf inf
14 c0_FFCA 0.248 0 0 0.248 False -inf inf
15 c0_FDCA 0.112 0 0 0.112 False -inf inf
16 c0_D_HMF 0 0 nan 0 False -inf inf
17 c0_D_DFF 0 0 nan 0 False -inf inf
18 c0_D_HMFCA 0 0 nan 0 False -inf inf
19 c0_D_FFCA 0 0 nan 0 False -inf inf
20 c0_D_FDCA 0 0 nan 0 False -inf inf
21 c0_Dx_HMF 0 0 nan 0 False -inf inf
22 c0_Dx_DFF 0 0 nan 0 False -inf inf
23 c0_Dx_HMFCA 0 0 nan 0 False -inf inf
24 c0_Dx_FFCA 0 0 nan 0 False -inf inf
25 c0_Dx_FDCA 0 0 nan 0 False -inf inf

To assess the fit results the evolution of concentration over time and charge passed over time are plotted, using respectively the functions plot.plot_c() and plot.plot_q(). For the concentrations evolution over time the lines correspond to the fit result and the points with errorbars to the experimental data.

from chemical_kinetics import plot

plot.plot_c(ds, ["HMF"])
_images/HMF_oxidation_WO3_37_0.svg
plot.plot_c(ds, ["DFF", "HMFCA", "FFCA", "FDCA"])
_images/HMF_oxidation_WO3_38_0.svg
plot.plot_q(ds)
_images/HMF_oxidation_WO3_39_0.svg