\documentclass[../book.tex]{subfiles}
\graphicspath{{\subfix{.}}}
\begin{document}

\chapter{Simulating Chemical Kinetics}
\label{07SimKin}

\setcounter{section}{-1}  % Start subsections from 0

\input{prereqs/07SimKin_prereqs}

\newpage
\section{Problem}
\setcounter{subsection}{-1}  % Start subsections from 0

With all the extra time you have gained automating your problem-solving with coding, you turn back to thinking about your project.  It is nice to be able to do these high throughput experiments, but how do you know what range of concentrations to start with, and how long the reaction will take to run to completion?  You would like to be able to predict what the concentrations of all of the components of the reaction are, not only to understand what is happening, but to be able to design experiments to probe the mechanism in more detail.

You think, ``the mechanism for this reaction is complicated, none of the solutions that I learned in p-chem class will work here.'' You also think about how useful it was to simulate titration behavior. ``Maybe I can write a program that simulates the reaction without having to do math that may be complicated or impossible?''

The overall reaction is simple, just \ce{A + B ->[$k_{cat}$] C}.  However, there are two suspected intermediates, \ce{X} and \ce{Y}, making the mechanism more complex:
\begin{center}
  
    \ce{A + cat ->[$k_1$] X}

    \ce{X + B ->[$k_2$] Y}

    \ce{Y ->[$k_3$] C + cat}
  
\end{center}


To figure out where to even start, you check a kinetic textbook out of the library.\footnote{A good reference is Steinfeld, Francisco, and and Hase's \textit{Chemical Kinetics and Dynamics}.  See a discussion about numerical solutions to kinetics problems in Chapter 2.6.} You find that chemical kinetics are an example of a initial value problem: you only know the concentrations of the species in your reaction at the beginning of the reaction.  However, if you can work out the differential rate laws, then you also know how the change of the concentrations vs. time depend on concentration.  So, you should be able to use these rate laws to project the behavior of the reaction forward in time.  The process of doing this is solving the set of differential equations, and luckily, chemical kinetics are usually represented by ordinary differential equations that can be solved straightforwardly (in principle).

You were taught to solve kinetics problems by first writing down differential rate laws for each component. For the reaction under consideration, you come up with:
\begin{equation}
\begin{split}
    \diff{\ce{[A]}}{t} & = -k_1 \ce{[A][cat]} \\
    \diff{\ce{[X]}}{t} & = k_1 \ce{[A][cat]} - k_2 \ce{[X][B]} \\
    \diff{\ce{[B]}}{t} & = -k_2 \ce{[X][B]} \\
    \diff{\ce{[Y]}}{t} & = k_2 \ce{[X][B]} - k_3 \ce{[Y]} \\
    \diff{\ce{[C]}}{t} & = k_3 \ce{[Y]} \\
    \diff{\ce{[cat]}}{t} & = -k_1 \ce{[A][cat]} + k_3 \ce{[Y]} \\
\end{split}
\end{equation}

You \textit{might} be able to make some approximations and solve this \cite{analytical}ly in some specific cases, but you also suspect that something in the solution might be inhibiting one of the steps of the cycle and would like to be able to account for that, too.  This means things will get even more complex, and it will probably be easier to just do this \cite{numerical}ly.  So, armed with your kinetics book and the internet, you set out to learn how to do it.

\section{Solution}
\setcounter{subsection}{-1}  % Start subsections from 0

A quick search shows that the Python \cite{library} \lstinline{scipy} can integrate systems of ordinary differential equations (ODEs) (among many other things!).  You see that there are different types of ODEs and a number of \cite{function}s to solve them built into the \lstinline{scipy.integrate} \cite{module}.  Kinetics are a type of initial value problem (you know the initial conditions only), rather than a boundary value problem (in which you know the value of the resulting function at both the beginning and end of the integration).  How do you map this onto a numerical solution?  

The rate laws are ODEs that tell you how fast the concentrations are changing in time as a function of all concentrations.  So, if you know all of the concentrations at any time, $t$, then you can make an estimate of the values of the concentrations at time $t + \Delta{t}$ by integrating the rates of change over that $\Delta{t}$.  You can then use the concentrations at $t + \Delta{t}$ to integrate up to $t + 2\Delta{t}$.  If you keep following the same process, you will have an estimate of the concentrations at a number of points in time in the future, from which you can make the familiar kinetics traces.  Smaller values of $\Delta{t}$ will tend to provide better estimates of the true solution, but require more computation.

Luckily for you, you don't have to figure out how to do the integration, just how to use an integration routine provided by \lstinline{scipy}!  %The main thing to know is that most ODE solvers choose the size of $\Delta{t}$ dynamically for each step of the process (known as adaptive solvers), which preserves accuracy while reducing the amount of computation (and thus, time) to reach the solution.

\subsection{Simulating the kinetics of one elementary reaction using \lstinline{scipy.solve_ivp()}}
\label{Simulating_the_kinetics_of_one_elementary_reaction_using_scipy.solve_ivp}

Looking through the documentation of the \lstinline{scipy.integrate} \cite{module}, it looks like the \cite{function} \lstinline{solve_ivp()} does what you want---determine the value(s) of one or a set of quantities as a function of time given initial values and ODEs for the rates of change.\footnote{For more information on numerical solutions of ordinary differential equations (particularly methods to evaluate accuracy when you don't already know the right answer) check out books like Hildebrand, F. B. \textit{Introduction to Numerical Analysis}; Dover Press, 1987, or Eric Sullivan's \textit{Numerical Methods: An Inquiry-Based Approach With Python} (\href{https://numericalmethodssullivan.github.io/}{https://numericalmethodssullivan.github.io}).}  However, you need to translate the discussion of the usage of the function to familiar chemistry terms.  It becomes clear that the $y$-value in the function will be your concentration, and \lstinline{y0} will be the initial concentration.  How to include the rate law is a little confusing, though.  \lstinline{solve_ivp()} has as the first required \cite{argument} \lstinline{fun}, which is a \cite{function} that contains the time derivative of $y$ and should be of the form \lstinline{fun(t, y)}.  The function \lstinline{fun} then needs to return the value of $\textrm{d}y/\textrm{d}t$.  It looks like you are required to have the \lstinline{t} argument there, even though your differential equations will not explicitly depend on time.  The second required argument is \lstinline{t_span}, which is a 2-element list or array that contains the start and stop time for the reaction simulation.  The last required argument is the initial value \lstinline{y0}.  You also see that there are a number of \cite{keyword arguments} that are optional.  Several appear to be notable:
\begin{itemize}
    \item \lstinline{method}: which algorithm to use to solve the IVP.  They are primarily broken down into methods better used on ``stiff" problems (of which many kinetics problems are examples) and ``non-stiff" problems (only very simple kinetics problems qualify).  The basic difference between these cases is that stiff problems have some components (reactants or intermediates) that may change very quickly with time and others that change more slowly.  Stiff solvers can solve these problems without resorting to extremely small step sizes, leading to more accurate results more quickly.  Typically, the default solver will be fine.
    
    \item \lstinline{args}: additional arguments to pass to your differential rate law function beyond the required \lstinline{t} and \lstinline{y}.  In your case, you don't want to hard code the rate constants into the rate laws function because you would have to change the values in the function every time.  Instead, you can make your function take an argument for \lstinline{k} as well and pass it whatever rates you want using this keyword.
    
    \item \lstinline{options}: parameters that mostly control the speed and accuracy of the solution.  Of these, the ones you are most likely to use are \lstinline{rtol} and \lstinline{atol}, which set the relative and absolute tolerances of the solution, respectively.  Lower values of both will yield more accurate solutions but require more time to solve.  \lstinline{rtol} basically controls the fractional accuracy, i.e. the percent difference, while \lstinline{atol} controls the absolute difference.  Thus, if your concentration is something like $1\times10^{-3}$ M, then \lstinline{rtol = 1e-3} will request accuracy to a part per thousand while \lstinline{atol = 1e-3} will request accuracy to the first decimal place.
\end{itemize}

Everything that you could want to know about the solution to the ODE, particularly the time points and the values of the function at each time, is \cite{return}ed in a single \cite{object}.  It contains several notable \cite{attribute}s:
\begin{itemize}
    \item \lstinline{t}: the time points at which the solution is known.
    
    \item \lstinline{y}: the values of the solution at each time point.  This is a 2-D \cite{Numpy_array}, where one dimension denotes the component of the reaction (each individual species present in the mechanism) and the second is corresponding concentration at each time point.
    
    \item \lstinline{status}: how the solution ended up, which can be used to make sure everything went okay.  This is 0 if the solver finished with no issues, and 1 or -1 in the rare case that something happened along the way that caused it not to finish.
\end{itemize}

\subsection{Planning}

Mapping this functions usage onto your problem, you decide that you will need the following:
\begin{itemize}
    \item Initial concentrations of all components of the reaction
    \item All of the rate constants
    \item A function containing the differential rate laws
    \item The total amount of time that you want to simulate the reaction
\end{itemize}
Thinking through what is required for a solution, you sketch out the four steps you need to solve the full problem:
\begin{enumerate}
    \item \textit{Define \cite{function}s containing the differential rate laws that need to be integrated.}

    The function should mimic the general form of the differential rate laws,
    \begin{equation}
        \diff{y}{t} = f(t, y, k),
    \end{equation}
    so it will need the \cite{argument}s \lstinline{t}, \lstinline{y}, and \lstinline{k} and \cite{return} a number $\diff{y}{t}$.  If you have more than one species in your reaction, this function will need to take an array of \lstinline{y}s and \lstinline{k}s and return an array of $\diff{y}{t}$s.  Because of this, even if you have only one component in your reaction, you will still need to provide a one-item list or Numpy array of \lstinline{y0} and will get back a one-item list with $\diff{y}{t}$.
    
    \item \textit{Set up the initial concentrations and the rates.}

    The initial concentrations will be stored as a list and passed as positional arguments to \lstinline{solve_ivp()}.  The rates will also be stored as a list and passed through \lstinline{solve_ivp()} to the differential rate laws function using the \lstinline{args} keyword.
    
    \item \textit{Invoke \lstinline{solve_ivp()}}

    This will just be a single line passing teh needed arguments to \lstinline{solve_ivp()} and assigning the results to a new variable that holds the object containing all the information about the solution.
    
    \item \textit{Examine and plot the results}

    You can then unpack only the information you need from the results object.  Specifically, you will want to plot all of the concentrations with respect to time!

\end{enumerate}
That doesn't seem too bad, but it is probably best to start with something simple where you already know the answer, just to make sure you understand the coding part.

\subsection{Validating coding approaches with simple test cases}
\label{Validating_coding_approaches_with_simple_test cases}


You figure that the easiest thing you can do is to model a simple unimolecular reaction,
\begin{center}
    \ce{A ->[$k$] B}
\end{center}
with the differential rate law $\diff{\ce{[A]}}{t} = -kt$.  You remember that  this type of reaction is described analytically by: $[\ce{A}](t) = [\ce{A}_0] \exp{(-kt)}$. Thus, you can use \lstinline{solve_ivp()} and  test the accuracy of the solutions you obtain from it against this analytical equation.

Now you can sketch out a plan for the solution:
\begin{figure}[H]
    \centering
    \includegraphics[]{ChapNarrative/07SimKin/Figures/FlowChartCh8.pdf}
    \caption{Flowchart for basic simulation of chemical kinetics using \lstinline{solve_ivp()} from scikit.}
    \label{fig:SimKin:ActivitDiagramBasicSim}
\end{figure}

\lstinputlisting[language=Python]{../../Distribution/Code by Chapter/Ch08 - Simulating Chemical Kinetics/FirstOrderKineticsSolver.py}

\subsection{Line-by-line explanation}

\begin{footnotesize}
\begin{itemize}
    \item [7--9:] We need to import \lstinline{numpy}, \lstinline{scipy.integrate} (which contains \lstinline{solve_ivp()}) and \lstinline{make_subplots} from \lstinline{plotly.subplots} to plot the results.
    
    \item [12--26:] This is the definition of the differential rate law. \lstinline{solve_ivp()} has the required parameter \lstinline{t}, even though it does not appear in the function.  The required parameter \lstinline{y} will hold the concentration value, and \lstinline{k} the rate.  We simply calculate the derivative of concentration with respect to time using the equation for the differential rate law and return its value (\lstinline{dydt}).

    \item [29--31:] We define the values of the initial concentration, the rate, and the beginning and end of the simulation time.  We don't have to do it on separate lines, we could just input their values as literals in the positional arguments and the \lstinline{args} keyword when we call \lstinline{solve_ivp()}.  However, it is clearer to do it separately and it makes it easier to reuse these values later in the code without having to remember to change them everywhere we use them.

    \item [34:] We call \lstinline{solve_ivp()} here and pass the necessary arguments and keywords to start solving the ODE.  Once execution finishes, \lstinline{solve_ivp()} returns an object containing arrays holding the \lstinline{t} and \lstinline{y} values for the solution and information to help diagnose the solution, which we store in the new variable \lstinline{solution}.

    \item [37:] Now we need to compute the analytical (correct) values for concentration using the known solution for first-order reactions at each time point in \lstinline{t} that we have numerical values.  The analytical concentrations are then stored in a new variable \lstinline{analytical_y}.  This will allow us to directly compute the absolute and relative errors of the numerical solution.

    \item [40--56:] Finally, we make a figure with two subplots.  The top subplot displays the analytical and numerical solutions for the concentration vs. time, and the bottom subplot displays the relative error (percent difference) between the two solutions.  
    
\end{itemize}
\end{footnotesize}

Running this code produces the following output:
\begin{figure}[H]
    \centering
    \includegraphics[width=4in]{ChapNarrative/07SimKin/Figures/FirstOrderExample.png}
    \caption{The result of your code to solve a simple first order reaction.  (top) A plot of concentration vs. time for the numerical solution and the known analytical one.  (bottom) The percent difference between the numerical and ``correct'' solutions.}
    \label{fig:SimKin:BasicResult}
\end{figure}

The first thing you notice is that the time points are not evenly spaced for the numerical solution.  Digging into the documentation for \lstinline{solve_ivp()}, you find that the default solver, ``RK45,'' is an adaptive step size solver, which means that it takes as big of step sizes as it can while trying to balance accuracy.  You also notice that the percent difference between the numerical solution and correct one increases with simulation time.  This is likely because the concentration at long times is close to zero, so the denominator of the percent difference becomes very small, blowing up the value.  Simply changing from percent difference to absolute difference shows this is the case:
\begin{figure}[H]
    \centering
    \includegraphics[width=4in]{ChapNarrative/07SimKin/Figures/FirstOrderExample-absolute.png}
    \caption{Like Figure \ref{fig:SimKin:BasicResult} but plotting the absolute rather than percent difference.}
    \label{fig:SimKin:ResultPercentage}
\end{figure}

The simulation is accurate to the order of 100s of nM or less, not bad!  You can use the keywords \lstinline{rtol} and \lstinline{atol} to increase the accuracy of these simulations if needed.

\subsection{Simulating multi-step reaction mechanisms using \lstinline{scipy.solve_ivp()}}
\label{Simulating_multi-step_reaction_mechanisms_using_scipy.solve_ivp}
With this example understood, you feel confident that you are ready to tackle your more complex mechanism.  You realize that you can just make small modifications to your code: you need to change the function that contains the differential rate equations, you need to have six initial concentrations and three rates, and you need to change the plotting code.  Other than that, you're all set.  

\section{Code}

You produce the following code:

\lstinputlisting[language=Python]{../../Distribution/Code by Chapter/Ch08 - Simulating Chemical Kinetics/CatalysisKineticsSolver.py}

Running this code gives you the following plot:
\begin{figure}[h]
    \centering
    \includegraphics[width=4in]{ChapNarrative/07SimKin/Figures/CatalysisExample.png}
    \caption{The output of the simulation using the code above to simulate the full kinetics of the chemical reactions outlined at the top of the chapter.}
    \label{fig:SimKin:FullSimPlot}
\end{figure}

You realize you can now solve any kinetics problem for which you can write the differential rate laws, and wonder why they didn't just teach you this in p-chem class... 

\subsection{Line-by-line explanation}

\begin{footnotesize}
\begin{itemize}
\footnotesize
    \item Most of the code is the same as the previous one.  Specific differences are:

    \item [27:] Since we still only have one variable, \lstinline{y}, we are passing the concentrations as a list.  For clarity, we are first unpacking the list into variables with helpful names.

    \item [28:] We are in the same situation with the rates, so we are unpacking them here, too.

    \item [30--35:] Now we are defining our system of rate laws instead of just one.

    \item [37:] We need to return the derivatives to \lstinline{solve_ivp()} as a list (or array) in the same order as the \lstinline{y}-values.

    \item [48--53:] Since we don't have analytical solutions, there's no way to plot the error.  Instead, we are just plotting all of the concentrations with respect to time using \lstinline{quickScatter()}.  Here we make use of built-in functionality in \lstinline{quickScatter()} whereby we can pass a list or array of $y$-arrays and it will automatically plot them each as a separate trace.  We also pass a list as the \lstinline{name} argument with names for each corresponding $y$-array.
    
\end{itemize}
\end{footnotesize}


% \section{Adding a Graphical User Interface}

% The main thing that is irritating you now is having to rerun your script every time you want to change a rate or initial concentration.  It would be nice for this to be more interactive.  You think about the web apps that you're used to using and wonder how hard it would be to make one of your own.  Looking into it, you realize that what you want is a \textbf{Graphical User Interface}, and that there's a bewildering array of options that you could look into.  You can check out some options and examples on the companion website for this book.

\hrule
\input{Exercises/Chapter 08 SimKin}
\hrule

\end{document}
