\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[english]{babel}



\usepackage{float}








\iflatexml
\documentclass[12pt]{article}
\fi

  \usepackage{fixltx2e}

\usepackage[top=1in, bottom=1.in, left=1in, right=1in]{geometry}

\usepackage{graphicx}

\usepackage{tabularx}

\usepackage{amssymb}

\usepackage{amsmath}

% \usepackage[version=3]{mhchem} % Not (yet) supported by LaTeXML.

\usepackage[numbers, super, sort&compress]{natbib}

\usepackage{setspace}

\usepackage[linktocpage,pdfstartview=FitH,colorlinks,linkcolor=blue,anchorcolor=blue,citecolor=blue,filecolor=blue,menucolor=blue,urlcolor=blue]{hyperref}


\begin{document}

\title{Evaluation of the Degree of Rate Control via Automatic Differentiation}



\author[1]{Yilin Yang}%
\author[1]{Siddarth Achar}%
\author[1]{John Kitchin}%
\affil[1]{Carnegie Mellon University}%


\vspace{-1em}



  
  \date{May 19, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
The degree of rate control quantitatively identifies the kinetically relevant (sometimes known as rate-limiting) steps of a complex reaction network. This concept relies on derivatives which are commonly implemented numerically, e.g. with finite differences. Numerical derivatives are tedious to implement, and can be problematic, and  unstable or unreliable. In this work, we demonstrate the use of automatic differentiation in the evaluation of the degree of rate control. Automatic differentiation libraries are increasingly available through modern machine learning frameworks. Compared to the finite differences, automatic differentiation provides solutions with higher accuracy with lower computational cost.  Furthermore, we illustrate a hybrid local-global sensitivity analysis method, the distributed evaluation of local sensitivity analysis (DELSA), to assess the importance of kinetic parameters over an uncertain space. This method also benefits from automatic differentiation to obtain high-quality results efficiently.%
\end{abstract}%



\sloppy


\doublespace

\section*{Introduction}
\label{sec:org0a451eb}

Complex reaction systems involving multiple reaction steps are common in the real world. Quantifying the influence of each reaction in a complex system is an important step to better understand what controls the kinetic behavior of the overall rate or to improve the reaction system to meet a desired target by adjusting the reaction conditions or catalysts. For example, we can increase the rate of the rate-determining step to increase the net reaction rate. Among various tools, the degree of rate control (DRC) is a versatile concept proposed by Campbell to measure the kinetic contribution of each reaction steps to a target reaction rate \hyperref[csl:1]{(Carsten Stegelmann \& Campbell, 8077}; \hyperref[csl:2]{Campbell, 2770)}.  There are many applications of this concept in the research areas of catalysis and microkinetic modeling \hyperref[csl:3]{(C.~Stegelmann \& P.~Stoltze, 2004}; \hyperref[csl:4]{\textit{{Mechanism of the Water Gas Shift Reaction on Pt : First Principles Experiments, and Microkinetic Modeling}}, 4608}; \hyperref[csl:5]{Christopher~A. Wolcott \& Campbell, 2015}; \hyperref[csl:6]{J{\o}rgensen \& Gr{\"o}nbeck, 4034}; \hyperref[csl:7]{Mao \& Campbell, 9465)}. These applications involve the investigation of the relationship between the reaction conditions and the rate-determining step \hyperref[csl:3]{(C.~Stegelmann \& P.~Stoltze, 2004)}, using the degree of rate control to screen the catalysts \hyperref[csl:5]{(Christopher~A. Wolcott \& Campbell, 2015)}, exploring the mechanism of reactions \hyperref[csl:3]{(C.~Stegelmann \& P.~Stoltze, 2004}; \hyperref[csl:4]{\textit{{Mechanism of the Water Gas Shift Reaction on Pt : First Principles Experiments, and Microkinetic Modeling}}, 4608)} and so forth. The kinetic DRC is defined as the derivative of the rate with respect to the energy of the transition state or the forward kinetic rate constant given the corresponding equilibrium constant is fixed:

\begin{equation}
\label{k-DRC} 
X_{RC, i} = \frac{k_i}{r} \left( \frac{\partial r}{\partial k_i} \right)_{k_{j \neq i}, K_i} = \left( \frac{\partial \ln{r}}{\partial \ln{k_i} } \right)_{k_{j \neq i}, K_i}

\end{equation}

where \(r\) is the net rate of the production of interest, \(k_i\) is the forward kinetic constant of step \(i\) and \(K_i\) is the equilibrium constant of step \(i\). The thermodynamic version of the DRC extends the application of this concept to the free energy of the intermediate species in the reaction system, which is mathematically defined as

\begin{equation}
X_{TRC, n}  = \left( \frac{\partial \ln{r}}{\partial \frac{-G_n^0}{RT} } \right)_{G_{m \neq n}^0, G_i^{0, TS}}

\end{equation}

where \(G^{0}_m\) is the free energy of species \(m\), \(G^{0, TS}_{i}\) is the transition-state energy of step \{i\}. Several variants of the DRC have been proposed to fulfill different purposes, including the DRC for selectivity \hyperref[csl:8]{(Talin Avanesian \& Christopher, 2016)}, for transient kinetics, \hyperref[csl:9]{(Campbell, 2001}; \hyperref[csl:2]{Campbell, 2770}; \hyperref[csl:10]{Foley \& Bhan, 2556)} and for uncertain parameters \hyperref[csl:11]{(Tian \& Rangarajan, 1353)}.

From a practical perspective, a simple and common way to assess the DRC is using finite difference (FD) approximations for the derivatives  \hyperref[csl:9]{(Campbell, 2001}; \hyperref[csl:10]{Foley \& Bhan, 2556}; \hyperref[csl:12]{\textit{{Quantifying the Impact of Parametric Uncertainty on Automatic Mechanism Generation for \ce Co2 Hydrogenation on Ni (111), Apr}}, 1437)}. These are fairly straightforward to implement and only require a few more lines of code in addition to the original simulation code. Mathematically, the centered difference approximation is formulated by

\begin{equation}
X_{RC, i}  =  \frac{k_i}{r} \left( \frac{\partial r}{\partial k_i} \right)_{k_{j \neq i}, K_i} \approx \frac{k_i}{r}\frac{\delta r}{\delta k_i}

\end{equation}

where \(\delta k_i\) is the perturbation applied on \(k_i\) and \(\delta r\) is the change of the net reaction rate resulted from the perturbation of \(k_i\). Although the finite difference approximation is a popular choice because of its ease of understanding and implementation, one must be careful to choose the magnitude of the perturbation. The change of the kinetic constant should be small enough such that the response of the net reaction is linear, but not so small that goes beyond the precision limitation of the computer \hyperref[csl:9]{(Campbell, 2001)}.

Formally, there is a trade-off between the truncation error and the rounding error in this approach. The scale of the truncation error is \(O(\delta)\) for the first-order derivative, which prefers small \(\delta\). However, when the \(\delta\) is smaller than the precision limitation, then the value is no longer reliable. In addition, FD requires \(O(n)\) rounds of function calls or forward simulations to get the derivatives of \(n\) parameters, which is time consuming for large size of parameters. For example, Bjarne et al. took 300 CPU-hours to conduct the sensitivity analysis when investigating the mechanism of \ce{CO_2} hydrogenation on Ni(111) using finite difference method \hyperref[csl:12]{(\textit{{Quantifying the Impact of Parametric Uncertainty on Automatic Mechanism Generation for \ce Co2 Hydrogenation on Ni (111), Apr}}, 1437)}.

To avoid these issues, sensitivity analysis methods like the direct sensitivity analysis and the adjoint sensitivity analysis \hyperref[csl:13]{(Yang Cao \& Serban, 1076)} are usually adopted by common differential equation solvers \hyperref[csl:14]{(\textit{{Sundials}}, 2005}; \hyperref[csl:15]{Zhang \& Sandu, 2014}; (missing citation) to provide the derivative of the numerical solution to the parameters for differential equations. The direct sensitivity analysis converts the solution sensitivity with respect to the parameters of differential equations into  \(n\) extra (number of parameters) differential equations, which are solved simultaneously with the original differential equations. The adjoint sensitivity analysis requires the definition of some scalar functional of the numerical solution and the parameters. Then the sensitivities are given by an integration. For the mathematical details, one could refer to the introduction of these two methods in the chemical kinetic systems by Sandu \hyperref[csl:17]{(Adrian Sandu \& Carmichael, 5083)}.

To be applied in the calculation of the DRC, the partial derivatives of the reaction rate to the concentrations and the kinetic parameters still needs to be solved, since the sensitivity analysis methods only provide the derivatives of the numerical solution (commonly the concentration or coverages in the chemical kinetic systems) to the parameters. One of the tools to integrate these derivatives is  automatic differentiation (AD), \hyperref[csl:18]{(Rall, 1981}; \hyperref[csl:19]{Atilim~Gunes Baydin \& Siskind, 2015)} which automatically evaluates the derivatives of a function that is built on a set of atomic operations and functions (e.g., addition, multiplication, exp, log, etc.). The derivatives are generated by chain rule based on the derivatives of these elementary operations. For the DRC case, the numerical integration and post functions from the kinetic constants to the reaction rate could be regarded as a sequence of the atomic operations and functions. Thus, the derivatives of the reaction rate to the kinetic constants could be evaluated by the chain rule applied on this sequence of basic operations. In addition, AD calculates the derivatives simultaneously with the function evaluations, which makes it more efficient compared to the finite difference or direct sensitivity analysis.

In this work, we adopt  AD to the evaluation of the DRC. In section \ref{sec-ad}, we introduce the working mechanism of the AD to obtain the derivatives of a function automatically. In section \ref{sec-results}, we illustrate three case studies to check the correctness of the AD and to show its advantages over FD. Specifically, we take the hypothetical reaction scheme from Foley's non-steady DRC work as a simple case \hyperref[csl:10]{(Foley \& Bhan, 2556)}. For a slightly more complicated example, we use the water-gas shift reaction adopted from Motagamwala's maximum rate work \hyperref[csl:20]{(Motagamwala \& Dumesic, 2879)}. Finally, we show the application of AD to calculate the DRC of a more complex mechanism, the propylene partial oxidation on Cu(100), \hyperref[csl:21]{(Song \& Wang, 2743)} which involves 17 elementary steps. Further more, we demonstrate the utilization of the distributed evaluation of local sensitivity analysis (DELSA) \hyperref[csl:22]{(\textit{{Distributed Evaluation of Local Sensitivity Analysis ( Delsa ), with Application to Hydrologic Models}}, 2014)} to deal with the uncertain range of the kinetic parameters, which is a hybrid local-global sensitivity analysis method to identify the important parameters and the importance distribution over an uncertain range.

\section*{Methodology}
\label{sec:orgf900477}

\subsection*{Automatic Differentiation \label{sec-ad}}
\label{sec:org7823da9}

Automatic differentiation has two modes to generate the derivatives: the forward mode and the reverse mode \hyperref[csl:18]{(Rall, 1981}; \hyperref[csl:19]{Atilim~Gunes Baydin \& Siskind, 2015)}. In the forward mode, the computational graph starts with the input variables, and grows along the elementary operations and functions applied on the input variables. During the forward expansion, the function evaluations and derivative calculations take place simultaneously. In the reverse mode, there are two rounds of evaluations. The first one is the forward evaluation of the function values starts from the input variables. The second round is the back-propagation of the derivatives from the function output to the input.

We use the example \(y = \ln{(3x_1 - 2x_2)} + x_1 x_2\) with \((x_1, x_2) = (1, 1)\) as a prototype example to illustrate the workflow of the AD. The computational graph is shown in Figure \ref{ad-graph}. The details of the forward AD and the reverse mode are shown in Table \ref{fwd-ad} and Table \ref{rvs-ad}. The comparison between the forward and the reverse mode is clear in Table \ref{fwd-ad} and Table \ref{rvs-ad}. In the forward mode, all derivatives of the intermediate and the final results with respect to a specified input variable (\(x_1\) in this case) are calculated in one forward propagation. However, in the reverse mode, the derivative of a specified scalar output with respect to all the intermediate and input variables are obtained in one backward propagation. Thus, the forward mode is suitable for the functions with fewer input dimensions, while the reverse mode is more efficient for the functions with fewer output dimensions.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/01-comp-graph/01-comp-graph}
\caption{{\label{ad-graph}Computational graph for the example \(y = \ln{(3x_1 - 2x_2)} + x_1 x_2\). Each blue node contains the corresponding variable and the operation that applied on its parent node(s). For instance, \(w_3\) is the result of the multiplication between a constant \(3\) and another node variable \(w_1\).%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{table}[htbp]
\caption{{\label{fwd-ad}Forward AD for the example \(y = \ln{(3x_1 - 2x_2)} + x_1 x_2\). Left side shows the forward function evaluations. Right side shows the derivative of \(\frac{\partial y}{\partial x_1}\). Head dot means \(\frac{\partial}{\partial x_1}\)}}
\centering
\begin{tabular}{ll}
\hline
Function evaluation & Derivative calculation (from top to bottom)\\
\hline
\(w_1 = x_1 = 1\) & \(\dot{w_1} = \dot{x_1} = 1\)\\
\(w_2 = x_2 = 1\) & \(\dot{w_2} = \dot{x_2} = 0\)\\
\(w_3 = 3 w_1 = 3\) & \(\dot{w_3} = 3\dot{w_1} = 3\)\\
\(w_4 = 2 w_2 = 2\) & \(\dot{w_4} = 2\dot{w_2} = 0\)\\
\(w_5 = w_3 - w_4 = 1\) & \(\dot{w_5} = \dot{w_3} - \dot{w_4} = 3\)\\
\(w_6 = \ln{w_5} = 0\) & \(\dot{w_6} = \frac{1}{w_5} \dot{w_5} = 3\)\\
\(w_7 = w_1 w_2 = 1\) & \(\dot{w_7} = w_1 \dot{w_2} + \dot{w_1} w_2 = 1\)\\
\(w_8 = w_6 + w_7 = 1\) & \(\dot{w_8} = \dot{w_6} + \dot{w_7} = 4\)\\
\(y = w_8 = 1\) & \(\dot{y} = \dot{w_8} = 4\)\\
\hline
\end{tabular}

\end{table}\selectlanguage{english}
\begin{table}[htbp]
\caption{{\label{rvs-ad}Reverse AD for the example \(y = \ln{(3x_1 - 2x_2)} + x_1 x_2\). Left side shows the forward function evaluations. Right side shows the back-propagation of the derivative of \(\frac{\partial y}{\partial x_1}\).}}
\centering
\begin{tabular}{ll}
\hline
Function evaluation & Derivative calculation (from bottom to top)\\
\hline
\(w_1 = x_1 = 1\) & \(\frac{\partial y}{\partial x_1} = \frac{\partial y}{\partial w_1} \frac{\partial w_1}{\partial x_1} = 4\)\\
\(w_2 = x_2 = 1\) & \(\frac{\partial y}{\partial x_2} = \frac{\partial y}{\partial w_2} \frac{\partial w_2}{\partial x_2} = -1\)\\
\(w_3 = 3 w_1 = 3\) & \(\frac{\partial y}{\partial w_1} = \frac{\partial y}{\partial w_1} + \frac{\partial y}{\partial w_3} \frac{\partial w_3}{\partial w_1} = 4\)\\
\(w_4 = 2 w_2 = 2\) & \(\frac{\partial y}{\partial w_2} = \frac{\partial y}{\partial w_2} + \frac{\partial y}{\partial w_4} \frac{\partial w_4}{\partial w_2} = -1\)\\
\(w_5 = w_3 - w_4 = 1\) & \(\frac{\partial y}{\partial w_3} = \frac{\partial y}{\partial w_5} \frac{\partial w_5}{\partial w_3} = 1\)\\
 & \(\frac{\partial y}{\partial w_4} = \frac{\partial y}{\partial w_5} \frac{\partial w_5}{\partial w_4} = -1\)\\
\(w_6 = \ln{w_5} = 0\) & \(\frac{\partial y}{\partial w_5} = \frac{\partial y}{\partial w_6} \frac{\partial w_6}{\partial w_5} = 1\)\\
\(w_7 = w_1 w_2 = 1\) & \(\frac{\partial y}{\partial w_1} = \frac{\partial y}{\partial w_7} \frac{\partial w_7}{\partial w_1} = 1\)\\
 & \(\frac{\partial y}{\partial w_2} = \frac{\partial y}{\partial w_7} \frac{\partial w_7}{\partial w_2} = 1\)\\
\(w_8 = w_6 + w_7 = 1\) & \(\frac{\partial y}{\partial w_6} = \frac{\partial y}{\partial w_8} \frac{\partial w_8}{\partial w_6} = 1\)\\
 & \(\frac{\partial y}{\partial w_7} = \frac{\partial y}{\partial w_8} \frac{\partial w_8}{\partial w_7} = 1\)\\
\(y = w_8 = 1\) & \(\frac{\partial y}{\partial w_8} = \frac{\partial y}{\partial y} \frac{\partial y}{\partial w_8} = 1\)\\
\hline
\end{tabular}

\end{table}

Similar to the example above,  AD could also be applied on the numerical integration of the ODE systems. Consider the initial-value problem

\begin{equation}
\dot{\boldsymbol{\theta}} = \boldsymbol{f}(\boldsymbol{\theta}, \boldsymbol{k}, \boldsymbol{K}, t), \; \boldsymbol{\theta}(\boldsymbol{k}, \boldsymbol{K}, t=0) = \boldsymbol{\theta}_0

\end{equation}

where \(\boldsymbol{\theta}\) is the state vector, \boldsymbol{f} is the state derivative vector, \(\boldsymbol{k}\) and \(\boldsymbol{K}\) are the parameters of this ODE system, and \boldsymbol{\theta}\_0 is the initial state vector. For the sake of simplicity, we assume the explicit forward Euler method is used to solve this ODE. Thus, the update equation is

\begin{equation}
\boldsymbol{\theta}_{n+1} = \boldsymbol{\theta}_{n} + h \boldsymbol{f}(\boldsymbol{\theta}_{n}, \boldsymbol{k}, \boldsymbol{K}, t_n)

\end{equation}

where \(h\) is the step size. Therefore, to get the derivative of \(\boldsymbol{\theta}_{n+1}\) with respect to \(\boldsymbol{k}\), we have

\begin{equation}
\label{ode-ad} 
\frac{d \boldsymbol{\theta}_{n+1}}{d \boldsymbol{k}} = \frac{d \boldsymbol{\theta}_{n}}{d \boldsymbol{k}} + h \frac{d \boldsymbol{f}(\boldsymbol{\theta}_{n}, \boldsymbol{k}, \boldsymbol{K}, t_n)}{d \boldsymbol{k}}

\end{equation}

In the forward AD, the derivative of the terms at the RHS of Eq \ref{ode-ad} with respect to \(\boldsymbol{k}\) is propagated to the LHS. In the reverse AD, the derivative of \(\boldsymbol{\theta}_{n+1}\) with respect to the LHS is back propagated to the variables in the RHS. In both modes, the derivatives of the ODE solution with respect to the parameters could be obtained automatically. Even if there are post operations applied on the ODE solution (e.g., conversion of the concentrations to reaction rates), the derivatives of the final results to the parameters can be calculated in an end-to-end way as long as these operations are in the same computational graph. We note that some ODE solvers may have the step size (\(h\)) that is dependent on the parameters, in this case, the derivative of the step size with respect to the parameters are enforced to be zero during the implementation \hyperref[csl:23]{(Eberhard \& Bischof, 1999)}.

In the past decade, various AD packages have been developed for applications in machine learning \hyperref[csl:19]{(Atilim~Gunes Baydin \& Siskind, 2015)}. Typical examples include Pytorch \hyperref[csl:24]{(\textit{{Pytorch: an Imperative Style, High-Performance Deep Learning Library}}, 2019)} and Jax \hyperref[csl:25]{(\textit{{Jax : Composable Transformations of P Ython+ N Um P Y Programs}}, 2018)} in Python, and ForwardDiff \hyperref[csl:26]{(Jarrett Revels \& Papamarkou, 2016)} in Julia. To integrate the AD into the solution of an ODE system, the ODE solver should be compatible to these AD packages and the operations and functions in the ODE solver should be included in the computational graph of the AD packages. There are several modules satisfying these requirements such as the torchdiffeq \hyperref[csl:27]{(Ricky T.~Q. Chen \& Duvenaud, 2018)} for Pytorch, the DifferentialEquations \hyperref[csl:28]{(Rackauckas \& Nie, 2017)} for ForwardDiff and the PyBaMM \hyperref[csl:29]{(\textit{{Python Battery Mathematical Modelling ( Pybamm}}, 2020)} for Jax. In this work, we use the \texttt{ForwardDiff} and \texttt{DifferentialEquations} in the Julia language since they provides various differentiable ODE solvers that could handle various non-stiff and stiff problems.

\subsection*{Distributed Evaluation of Local Sensitivity Analysis}
\label{sec:org202449b}

Distributed evaluation of local sensitivity analysis (DELSA) \hyperref[csl:22]{(\textit{{Distributed Evaluation of Local Sensitivity Analysis ( Delsa ), with Application to Hydrologic Models}}, 2014)} is a hybrid local-global sensitivity analysis method to measure the distribution of parameter sensitivity across the parameter space with low computational cost. Basically, it is an extension to local sensitivity analysis that takes the uncertainty of the parameters into account. The importance of a parameter over the parameter space is measured by a local sensitivity statistic like the median or the mean of a set of samples drawn from the parameter space. Local sensitivity analysis is then conducted on each parameter sample. This makes DELSA much cheaper to get the detailed distribution of the importance over the parameter space than the Sobol's indices which is a popular global sensitivity analysis method based on the variance decomposition \hyperref[csl:30]{(Sobol′, 2001}; \hyperref[csl:31]{Saltelli, 2002)}. Mathematically, the first-order sensitivity measure for \(j^{th}\) parameter at sample \(i\) is defined as

\begin{equation}
\label{delsa-sensi} 
S_{ij} = \frac{\left| \frac{y_i}{\theta_{ij}} \right|^2 s_j^2}{V(y_i)}

\end{equation}

where \(S_{ij}\) is the sensitivity measure for the parameter \(\theta_{ij}\) at sample \(i\), \(y_i\) is the model output or an element of the model output of sample \(i\), \(s_j\) is the prior variance of \(\theta_j\), and \(V(y_i)\) is the total variance of the model output \(y_i\), which can be evaluated using the first-order-second-moment method \hyperref[csl:32]{(Hill \& Tiedeman, 2005)}:

\begin{equation}
V{(y_i)} = \left( \frac{\partial y_i}{\partial \boldsymbol{\theta}} \right)^T \left( \boldsymbol{X}^T \boldsymbol{\omega} \boldsymbol{X}  \right) \left( \frac{\partial y_i}{\partial \boldsymbol{\theta}} \right)

\end{equation}

where \(\boldsymbol{X}\) is a matrix of (\(n_{obs}\) + \(n_{prior}\)) rows and \(n_{param}\) columns. \(n_{obs}\), \(n_{prior}\), \(n_{param}\) are the number of observations, the number of prior information equations and the number of parameters respectively. In the application of DELSA, there is no observation and the prior information is the variance of each parameter. Therefore, \(n_{obs} \; = \; 0\) and \(n_{prior} \; = \; n_{param}\), and each row of \(\boldsymbol{X}\) has zeros except for one which indicates the parameter associated with the prior information. \(\boldsymbol{\omega}\) is a diagonal matrix contains the reciprocal of the prior variance of each parameter. More details of the structure of \(\boldsymbol{X}\) and \(\boldsymbol{\omega}\) could refer to the related works \hyperref[csl:32]{(Hill \& Tiedeman, 2005}; \hyperref[csl:33]{Petropoulos \& Srivastava, 2016}; \hyperref[csl:22]{\textit{{Distributed Evaluation of Local Sensitivity Analysis ( Delsa ), with Application to Hydrologic Models}}, 2014)}.

Intuitively, the sensitivity measure \(S_{ij}\) captures the contribution of \(\theta_{ij}\) to the total uncertainty of the output \(y_{i}\). The parameters with more contribution to the uncertainty are considered more important. This sensitivity measure gives similar importance evaluation to the Sobol's method \hyperref[csl:30]{(Sobol′, 2001)} for the uncertain parameters in the previous reports \hyperref[csl:34]{(Sobol' \& S.~Kucherenko, 1212}; \hyperref[csl:22]{\textit{{Distributed Evaluation of Local Sensitivity Analysis ( Delsa ), with Application to Hydrologic Models}}, 2014)}. In addition, DELSA could be easily integrated to the original local sensitivity implementation which is the DRC with automatic differentiation in our work.

\section*{Results and Discussion \label{sec-results}}
\label{sec:org410b5cd}

We firstly show that AD can reproduce the DRC results of previous reports \hyperref[csl:20]{(Motagamwala \& Dumesic, 2879}; \hyperref[csl:10]{Foley \& Bhan, 2556)} in the first two simple cases. Then, we compare the performance of the FD and the AD on another more complicated reaction mechanism of propylene oxidation \hyperref[csl:21]{(Song \& Wang, 2743)}. Lastly, we show the application of AD to deal with uncertain parameters by manually introducing an uncertain range for the kinetic constants in the propylene oxidation case.

\subsection*{Case I: Hypothetical Two-Step Reaction}
\label{sec:org11fb286}

We first consider the hypothetical reaction mechanism in Table \ref{case1-scheme} \hyperref[csl:10]{(Foley \& Bhan, 2556)}.\selectlanguage{english}
\begin{table}[htbp]
\caption{{\label{case1-scheme}Hypothetical Two-Step Catalytic Reaction (case I). \hyperref[csl:10]{(Foley \& Bhan, 2556)}}}
\centering
\begin{tabular}{rlll}
Step id & Elementary Step & \(k_i\) & \(k_{-i}\)\\
\hline
1 & \(A + * \leftrightarrow A*\) & \(10^{-5}\) & 0\\
2 & \(A* + B \rightarrow C + *\) & 1 & NA\\
\hline
\end{tabular}

\end{table}

which leads to the following differential equations (Eqs \ref{case1-diff1} - \ref{case1-diff2}):

\begin{equation}
\label{case1-diff1} 
\frac{1}{L}\frac{\theta_{A*}(t)}{dt} = k_1 a_A(t) \theta_{*}(t) - k_{-1} \theta_{A*}(t) - k_2 a_B(t) \theta_{A*}(t)

\end{equation}

\begin{equation}
\label{case1-diff2} 
\frac{1}{L}\frac{\theta_{*}(t)}{dt} = -k_1 a_A(t) \theta_{*}(t) + k_{-1} \theta_{A*}(t) + k_2 a_B(t) \theta_{A*}(t)

\end{equation}

where \(\theta_i\) are the coverages, \(a_i\) are the thermodynamic activities and \(L\) is the number of active sites. The net reaction rate is defined as the rate to produce \(C\) per active site:

\begin{equation}
\frac{r_C(t)}{L} = k_2 a_B(t) \theta_{A*}(t)

\end{equation}

The analytical solution for \(r_C(t)/L\) is

\begin{equation}
\frac{r_C(t)}{L} = \frac{k_1 a_A k_2 a_B}{k_1 a_A + k_{-1} + k_2 a_B}\left( 1 - e^{-(k_1 a_A + k_{-1} + k_2 a_B)t} \right) + k_2 a_B \theta_{A*, 0} e^{-(k_1 a_A + k_{-1} + k_2 a_B)t}

\end{equation}

where \(\theta_{A*, 0}\) is the coverage of \(A*\) at \(t \; = \; 0\).

According to the settings in the original paper \hyperref[csl:10]{(Foley \& Bhan, 2556)}, \(a_A \; = \; 1\), \(a_B(t < 0) \; = \; 1\) and \(a_B(t \geq 0) \; = \; 3\). Since \(k_{-1}\) is set as 0 and \(k_{1} \ll k_{2}\), thus

\begin{equation}
\frac{r_C(t)}{L} = (k_1 a_A)\left(1 - e^{-(k_2 a_B)t} \right) + k_2 a_B \theta_{A*, 0} e^{-(k_2 a_B)t}

\end{equation}
which has the corresponding DRC as:

\begin{equation}
\label{case1-drc1} 
X_{RC, 1} = \frac{k_1 L}{r_C(t)} a_A\left( 1-e^{-(k_2 a_B)t} \right)

\end{equation}

\begin{equation}
\label{case1-drc2} 
X_{RC, 2} = \frac{k_2 L}{r_C(t)} \left[k_1 a_A a_B t e^{-(k_2 a_B)t}   + a_B \theta_{A*, 0} e^{-(k_2 a_B)t} -k_2 a_B^2 \theta_{A*, 0}t e^{-(k_2 a_B)t}   \right]

\end{equation}

We note that this transient DRC is based on the definition of Eq \ref{k-DRC}, which is different from the modified version proposed by Bhan \hyperref[csl:10]{(Foley \& Bhan, 2556)}. Upon Eq \ref{case1-drc1} - \ref{case1-drc2}, notice that at \(t \; = \; 0^{+}\), \(X_{RC, 1}\) is \(0\) and \(X_{RC, 2}\) is \(1\). As \(t \rightarrow \infty\), \(X_{RC, 1}\) grows to \(1\) and \(X_{RC, 2}\) decreases to \(0\). The solutions from the analytical form and the AD are shown in Figure \ref{case1-drc-fig}, in which the solution of the AD perfectly matches the analytic solution over the whole time range. The DRC in the Figure \ref{case1-drc-fig} does not obey the sum of kinetic DRC equaling one because the rate during the transient process also depends on the time, apart from the kinetic constants. This issue is discussed in more detail in the non-steady DRC paper \hyperref[csl:10]{(Foley \& Bhan, 2556)}. The main message conveyed by this simple example is that the AD is practically equivalent to having the analytical derivatives to evaluate the transient DRC defined by Eq \ref{k-DRC}, but without the need to analytically derive the expressions or to approximate them with finite differences.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/02-case1-DRC/02-case1-DRC}
\caption{{\label{case1-drc-fig}Degree of rate control for the hypothetical reaction of case 1 evaluated by analytic solution and automatic differentiation.%
}}
\end{center}
\end{figure}

\subsection*{Case II: Redox mechanism for water-gas shift}
\label{sec:org9123628}

Case I was a hypothetical example and it only contained two steps with manually set kinetic parameters, which is relatively simple. In case II, we consider a more complicated reaction mechanism for the water-gas shift reaction \hyperref[csl:20]{(Motagamwala \& Dumesic, 2879)} which is listed in Table \ref{case2-scheme}. The pressures for the gas-phase species \(P_{CO}\), \(P_{H_2O}\), \(P_{H_2}\), \(P_{CO_2}\) are \(0.07\) atm, \(0.21\) atm, \(0.38\) atm, and \(0.085\) atm, respectively. The target net rate is of the production of the hydrogen, which could be expressed as (Eq \ref{case2-net-rate}):

\begin{equation}
\label{case2-net-rate} 
r_{H_2} = k_7 \theta_{H*}^2 - k_{-7} \frac{P_{H_2}}{P^{\standardstate}} \theta_{*}^2

\end{equation}

For this reaction system, the DRC evaluated by the AD are shown in Figure \ref{case2-drc-fig}. Among the 7 steps of this reaction mechanism, only two steps are identified as important to the net rate. Step 4 (dissociation of \(OH*\)) has a DRC of 0.88 while step 5 (formation of \(CO_2 *\)) has a DRC of 0.12. The other 5 steps have little influence on the net rate of the whole reaction system. This result is consistent to the calculation of the original paper \hyperref[csl:20]{(Motagamwala \& Dumesic, 2879)} and illustrates the reliability of the AD to evaluation the DRC of a moderately complicated reaction system.\selectlanguage{english}
\begin{table}[htbp]
\caption{{\label{case2-scheme}Redox mechanism for water-gas shift (case II). \hyperref[csl:20]{(Motagamwala \& Dumesic, 2879)}}}
\centering
\begin{tabular}{rlll}
Step id & Elementary Step & \(K_{eq, i}\) & \(k_{i}\)\\
\hline
1 & \(CO + * \leftrightarrow CO*\) & \(2.15 \times 10^2\) & \(1.33 \times 10^8\)\\
2 & \(H_2O + * \leftrightarrow H_2O*\) & \(5.93 \times 10^{-5}\) & \(2.01 \times 10^{11}\)\\
3 & \(H_2 O* + * \leftrightarrow H* + OH*\) & \(6.28 \times 10^{-2}\) & \(2.64 \times 10^6\)\\
4 & \(OH* + * \leftrightarrow H* + O*\) & \(1.18 \times 10^{-5}\) & \(5.24 \times 10^1\)\\
5 & \(CO* + O* \leftrightarrow CO_2 * + *\) & \(1.03 \times 10^3\) & \(2.05 \times 10^5\)\\
6 & \(CO_2 * \leftrightarrow CO_2 + *\) & \(1.92 \times 10^5\) & \(1.48 \times 10^{12}\)\\
7 & \(2H* \leftrightarrow H_2 + 2*\) & \(4.50 \times 10^1\) & \(5.32 \times 10^2\)\\
\hline
\end{tabular}

\end{table}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/03-case2-DRC/03-case2-DRC}
\caption{{\label{case2-drc-fig}Degree of rate control for the water-gas shift reaction (case II). \protect\hyperref[csl:20]{(Motagamwala \& Dumesic, 2879)}%
}}
\end{center}
\end{figure}

\subsection*{Case III: Propylene Partial Oxidation}
\label{sec:org23f4f49}

Our third case is the propylene partial oxidation \hyperref[csl:21]{(Song \& Wang, 2743)}, which could be an environmentally friendly route to produce propylene oxide. Its 17 elementary steps and corresponding kinetic constants are shown in Table \ref{case3-scheme}. The partial pressures for the gas-phase species \(P_{C_3H_6}\) and \(P_{O_2}\) are 0.1 bar and 0.05 bar respectively. The reaction simulation temperature is 500 K. The net rate of interest is the desorption of propylene oxide, which is the sum of the step 13 and 14 (Eq \ref{case3-net-rate}):

\begin{equation}
\label{case3-net-rate} 
r_{PO} = k_{13} \theta_{PO_1*} - k_{-13} \frac{P_{PO}}{P^{\standardstate}} \theta_{v*} + k_{14} \theta_{PO_2*} - k_{-14} \frac{P_{PO}}{P^{\standardstate}} \theta_{v*}

\end{equation}

The DRC results for this case are shown in Figure \ref{case3-drc-fig}, where the adsorption of propylene as the III type and the desorption of \(PO_2*\) have the major positive contribution to the generation of \(PO(g)\). These positive DRCs are reasonable since the adsorption of propylene provides the material to produce \(PO(g)\) and the desorption of \(PO_2*\) directly generates \(PO(g)\). On the other side, the desorption of \(O_2*\) hinders the production of \(PO(g)\) since this step results in more \(v*\) produced and more \(O*\) and \(O1*\) consumed, which benefits the negative direction of step 13 and 14.\selectlanguage{english}
\begin{table}[htbp]
\caption{{\label{case3-scheme}Elementary steps for the propylene partial oxidation (case III) \hyperref[csl:21]{(Song \& Wang, 2743)}. \(E_a\) and \(E_a^{-r}\) are the activation energies for the forward and reverse directions. \(A\) and \(A^{-r}\) are the pre-expotential factor of the forward and reverse reactions. \(v\) represents the oxygen vacancy.}}
\centering
\begin{tabular}{rlllll}
Step id & Elementary Step & \(E_a (eV)\) & \(A\) & \(E_a^{-r} (eV)\) & \(A^{-r  }\)\\
\hline
1 & \(C_3H_6(g) + * \leftrightarrow C_3H_6(I)*\) & \(0.00\) & \(1.87 \times 10^8\) & \(0.42\) & \(1.00 \times 10^{13}\)\\
2 & \(C_3H_6(g) + * \leftrightarrow C_3H_6(II)*\) & \(0.00\) & \(1.87 \times 10^8\) & \(0.57\) & \(1.00 \times 10^{13}\)\\
3 & \(C_3H_6(g) + * \leftrightarrow C_3H_6(III)*\) & \(0.00\) & \(1.87 \times 10^8\) & \(0.58\) & \(1.00 \times 10^{13}\)\\
4 & \(C_3H_6(I)* + O* \leftrightarrow C_3H_5* + OH*\) & \(0.36\) & \(8.03 \times 10^{12}\) & \(1.40\) & \(8.93 \times 10^{12}\)\\
5 & \(C_3H_6(II)* + O* \leftrightarrow OMP_1*\) & \(0.59\) & \(1.17 \times 10^{13}\) & \(1.10\) & \(1.78 \times 10^{13}\)\\
6 & \(C_3H_6(III) + O* \leftrightarrow OMP_2*\) & \(0.31\) & \(1.28 \times 10^{13}\) & \(0.95\) & \(1.92 \times 10^{13}\)\\
7 & \(C_3H_5* +O* \leftrightarrow C_3H_4O* + H_2O*\) & \(0.30\) & \(1.14 \times 10^{13}\) & \(1.19\) & \(7.14 \times 10^{12}\)\\
8 & \(C_3H_5O* + OH* \leftrightarrow C_3H_4O* + H_2O*\) & \(0.54\) & \(1.09 \times 10^{13}\) & \(1.72\) & \(1.54 \times 10^{13}\)\\
9 & \(OMP_1* \leftrightarrow PO_1* +*\) & \(0.62\) & \(4.10 \times 10^{13}\) & \(0.83\) & \(1.53 \times 10^{13}\)\\
10 & \(OMP_2* \leftrightarrow PO_2* + *\) & \(0.77\) & \(3.24 \times 10^{13}\) & \(0.89\) & \(1.31 \times 10^{13}\)\\
11 & \(H_2O* \leftrightarrow H_2O(g) + v*\) & \(0.76\) & \(1.00 \times 10^{13}\) & \(0.00\) & \(2.85 \times 10^8\)\\
12 & \(C_3H_4O* \leftrightarrow C_3H_4O(g) + v*\) & \(0.10\) & \(1.00 \times 10^{13}\) & \(0.00\) & \(1.62 \times 10^8\)\\
13 & \(PO_1* \leftrightarrow PO(g) + v*\) & \(0.90\) & \(1.00 \times 10^{13}\) & \(0.00\) & \(1.59 \times 10^8\)\\
14 & \(PO_2* \leftrightarrow PO(g) + v*\) & \(0.96\) & \(1.00 \times 10^{13}\) & \(0.00\) & \(1.59 \times 10^8\)\\
15 & \(O_2* + v* \leftrightarrow O* + O1*\) & \(0.00\) & \(1.00 \times 10^{13}\) & \(1.43\) & \(1.00 \times 10^{13}\)\\
16 & \(2O1* \leftrightarrow O_2* + *\) & \(0.00\) & \(1.00 \times 10^{13}\) & \(1.18\) & \(1.00 \times 10^{13}\)\\
17 & \(O_2* \leftrightarrow O_2(g) + *\) & \(1.36\) & \(1.00 \times 10^{13}\) & \(0.00\) & \(2.14 \times 10^8\)\\
\hline
\end{tabular}

\end{table}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/04-case3-DRC/04-case3-DRC}
\caption{{\label{case3-drc-fig}Degree of rate control for the propylene partial oxidation reaction (case III). \protect\hyperref[csl:21]{(Song \& Wang, 2743)}%
}}
\end{center}
\end{figure}

We then compare the DRC calculation results using FD and AD in this case. The temperature is set as 350 K in this comparison. We evaluate the transient DRC defined by Eq \ref{k-DRC} using the FD and the AD methods separately. The transient DRC for step 3, 4, 14, and 17 are shown in Figure \ref{case3-fd-ad}. Different perturbation magnitudes (\(10^{-4}\), \(10^{-11}\), \(10^{-14}\)) are used in the FD method. There is an optimal choice of the perturbation size to reduce the truncation error and the rounding error. A large perturbation size suffers from the truncation error resulting from the nonlinearity of the target function, which is the transient part from \(t \; = \; 0\) to \(t \; = \; 5\). Too small of a perturbation size makes the solution affected by the rounding error due to the limited precision for floating numbers. In this case, \(10^{-11}\) is the best perturbation size among these candidates, whose result is the most aligned with the solution of the AD method. Case III illustrates that although with finite difference it is possible to get a reliable transient DRC for a complicated system, it highly depends on the choice of the perturbation size. This issue does not hold for the automatic differentiation since there is no truncation error during the derivative evaluation process of the AD.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/05-case3-transient-350K/05-case3-transient-350K}
\caption{{\label{case3-fd-ad}Transient degree of rate control evaluated by the finite difference and automatic differentiation for the propylene partial oxidation reactions 3, 4, 14 and 17 (case III) \protect\hyperref[csl:21]{(Song \& Wang, 2743)}. Different perturbation sizes are used in the FD method as indicated in the legend.%
}}
\end{center}
\end{figure}

\subsection*{Degree of Rate Control for Uncertain Parameters}
\label{sec:orgaf5c742}

It is common to consider an uncertain space of the kinetic parameters instead of the exact values in the real-world catalyst applications. To illustrate the usage of the DELSA \hyperref[csl:22]{(\textit{{Distributed Evaluation of Local Sensitivity Analysis ( Delsa ), with Application to Hydrologic Models}}, 2014)} to measure the importance of the parameters over a space, we hypothetically add an uncertainty range to the forward kinetic constant of each elementary step of the case III, which corresponds to an uncertain range of [-0.03 eV, 0.03 eV] on the forward activation energy. We assume the log of the kinetic parameters are uniformly distributed in the candidate space. The samples are drawn using the quasi-random Sobol sequence \hyperref[csl:35]{(Sobol', 1967)} which more evenly samples the space than a uniform distribution would. According to the definition of the DELSA importance measure (Eq \ref{delsa-sensi}), the importance of a kinetic parameter \(k_{ij}\) to the net rate of sample \(i\) could be measured by:

\begin{equation}
S_{k_{ij}} = \frac{\left| \frac{\partial \ln{r_i}}{\partial \ln{k_{ij}}} \right|^2_{k_{im \neq ij}, K_{j}} s_{j}^2 }{V({\ln{r_i}})}

\end{equation}

where \(s_j\) is the variance of the uniform distribution of \(k_{j}\). Among the 17 parameters, only 3 parameters are found to be important (with the average DELSA importance measure \textgreater{} 0.1) over the potential space. Their distributions are shown in Figure \ref{case3-delsa}. It is notable that the important steps identified by DELSA (step 3, 14 and 17) are different from the local DRC (step 1, 4, 5, 17). This is because DELSA considers the contribution of a parameter over the whole space instead of a local point. An important parameter at local scope (e.g., step 1) could have small contribution from a global perspective. A trivial parameter locally (e.g., step 14) is possible to be a significant step globally. In addition, DELSA provides the distribution of the importance across the parameter space. For example, the histogram of step 14 has two peaks near 0 and 1, which means that for about half of the samples, step 14 has no effect on the net rate (with 0 importance measure), but for the other half of the samples, step 14 is the most important parameter (with 1 importance measure). High quality distributions rely on more samples in the potential space, which could be expensive for finite differences and difficult to ensure the results are accurate. Using the automatic differentiation, DELSA could be performed with higher accuracy and lower computational cost.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/06-case3-delsa/06-case3-delsa}
\caption{{\label{case3-delsa}DELSA results for 3 important steps (3, 14 and 17) in the propylene oxidation reaction \protect\hyperref[csl:21]{(Song \& Wang, 2743)}.%
}}
\end{center}
\end{figure}

\section*{Conclusion}
\label{sec:org7121d36}

The degree of rate control is a versatile concept in heterogeneous catalysis. We discussed the application of automatic differentiation (AD) in the evaluation of the derivatives needed to evaluate the DRC. We also compared the AD method to the commonly used finite difference (FD) method. Compared to FD, AD provides a faster and more accurate solution. There is no need to choose the optimal perturbation size in AD which is a critical step in FD to obtain a reliable result. In three cases with various complexity, we demonstrated the correctness of the AD method. In addition, we illustrated that AD could be used to perform the DELSA method, which is a hybrid local-global sensitivity analysis method to measure the importance of parameters over an uncertain space. With AD, DELSA could be conducted accurately and efficiently.

\section*{Acknowledgments}
\label{sec:org5fd39bb}

This material is based upon work supported under NSF DMREF Award CBET-1921946.

\newpage

\selectlanguage{english}
\FloatBarrier
\section*{References}\sloppy
\phantomsection
\label{csl:1}\textit{{Degree of Rate Control: How Much the Energies of Intermediates and Transition States Control Rates}}. (8077).

\phantomsection
\label{csl:2}\textit{{The Degree of Rate Control: A Powerful Tool for Catalysis Research}}. (2770).

\phantomsection
\label{csl:3}\textit{{Microkinetic Modeling of Ethylene Oxidation Over Silver}}. (2004).

\phantomsection
\label{csl:4}(4608).

\phantomsection
\label{csl:5}\textit{{Degree of Rate Control Approach to Computational Catalyst Screening}}. (2015).

\phantomsection
\label{csl:6}\textit{{Connection Between Macroscopic Kinetic Measurables and the Degree of Rate Control}}. (4034).

\phantomsection
\label{csl:7}\textit{{Apparent Activation Energies in Complex Reaction Mechanisms: A Simple Relationship Via Degrees of Rate Control}}. (9465).

\phantomsection
\label{csl:8}\textit{{Mechanism of \ce Co2 Reduction by \ce H2 on Ru and General Selectivity Descriptors for Late-Transition Metal Catalysts}}. (2016).

\phantomsection
\label{csl:9}\textit{{Finding the Rate-Determining Step in a Mechanism}}. (2001).

\phantomsection
\label{csl:10}\textit{{Degrees of Rate Control At Non(pseudo)steady-State Conditions}}. (2556).

\phantomsection
\label{csl:11}\textit{{Computing a Global Degree of Rate Control for Catalytic Systems}}. (1353).

\phantomsection
\label{csl:12}(1437).

\phantomsection
\label{csl:13}\textit{{Adjoint Sensitivity Analysis for Differential-Algebraic Equations the Adjoint Dae System and Its Numerical Solution}}. (1076).

\phantomsection
\label{csl:14}(2005).

\phantomsection
\label{csl:15}\textit{{Fatode : a Library for Forward, Adjoint, and Tangent Linear Integration of Odes}}. (2014).

\phantomsection
\label{csl:17}\textit{{Direct and Adjoint Sensitivity Analysis of Chemical Kinetic Systems with Kpp : Part I-Theory and Software Tools}}. (5083).

\phantomsection
\label{csl:18}\textit{{Emph Automatic Differentiation: Techniques and Applications}}. (1981).

\phantomsection
\label{csl:19}\textit{{Automatic Differentiation in Machine Learning: a Survey}}. (2015).

\phantomsection
\label{csl:20}\textit{{Analysis of Reaction Schemes Using Maximum Rates of Constituent Steps}}. (2879).

\phantomsection
\label{csl:21}\textit{{A Dft Study and Microkinetic Simulation of Propylene Partial Oxidation on Cuo (111) and Cuo (100) Surfaces}}. (2743).

\phantomsection
\label{csl:22}(2014).

\phantomsection
\label{csl:23}\textit{{Automatic Differentiation of Numerical Integration Algorithms}}. (1999).

\phantomsection
\label{csl:24}(2019).

\phantomsection
\label{csl:25}(2018).

\phantomsection
\label{csl:26}\textit{{Forward-Mode Automatic Differentiation in Julia}}. (2016).

\phantomsection
\label{csl:27}\textit{{Neural Ordinary Differential Equations}}. (2018).

\phantomsection
\label{csl:28}\textit{{Differentialequations.jl --A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia}}. (2017).

\phantomsection
\label{csl:29}(2020).

\phantomsection
\label{csl:30}\textit{{Global Sensitivity Indices for Nonlinear Mathematical Models and Their Monte Carlo Estimates}}. (2001).

\phantomsection
\label{csl:31}\textit{{Making Best Use of Model Evaluations to Compute Sensitivity Indices}}. (2002).

\phantomsection
\label{csl:32}\textit{{Emph Effective Groundwater Model Calibration: With Analysis of Data, Sensitivities, Predictions, and Uncertainty}}. (2005).

\phantomsection
\label{csl:33}\textit{{Emph Sensitivity Analysis in Earth Observation Modelling}}. (2016).

\phantomsection
\label{csl:34}\textit{{A New Derivative Based Importance Criterion for Groups of Variables and Its Link with the Global Sensitivity Indices}}. (1212).

\phantomsection
\label{csl:35}\textit{{On the Distribution of Points in a Cube and the Approximate Evaluation of Integrals}}. (1967).

\end{document}

