### Kinetic Models Describe the Folding Reaction
Our computational goal is to determine the underlying kinetic model by fitting time-resolved \(\cdot\)OH data. This data, however, does not directly measure the state curves \(\vec{x}(t)\). Instead, it measures the change in solvent accessibility of specific nucleotides during the folding reaction [14, 18, 19]. As the molecule folds, a subset of the nucleotides become progressively more buried (or more exposed) and therefore less reactive (more reactive) to the \(\cdot\)OH probe. We monitor the change in accessibility of these nucleotides as a function of time, which yields experimental time-progress curves (Figure 1C), which we call \(\vec{C}_{E}(t)\). As an example, we use data collected for the folding of the L-21 _T. thermophila_ group I intron [8] in the presence of 10 mM MgCl\({}_{2}\). The \(\cdot\)OH footprinting curves shown in Figure 1C correspond to individual subdomains of the molecule; in this case the green curve (Figure 1C) is the average change in accessibility of nucleotides in the P4P6 subdomain (Figure 1A, secondary structure), while the red curves correspond to the peripheral helices, and the blue curves correspond to nucleotides in the catalytic core [9]. These assignments are based on a k-means clustering of individual time-progress curves using the Gap statistic described previously [9]. The state curves, \(\vec{x}(t)\) (Figure 1B), and predicted \(\cdot\)OH footprinting, \(\vec{C}_{P}(t)\) (Figure 1C), curves are related by  
\[\vec{C}_{P}(t)=\mathcal{P}\vec{x}(t) \tag{6}\]  
where the matrix \(\mathcal{P}\) represents the linear combinations of the state curves (\(\vec{x}(t)\)) to generate the progress curves (\(\vec{C}_{P}(t)\)). This matrix has the dimensions of the number of states (\(U(t),I_{1}(t),\ldots,F(t)\)) by the number of number of intermediates plus one, which is also the length of \(\vec{C}_{P}(t)\). The matrix \(\mathcal{P}\) for the data set shown in Figure 1 is given by  
\[\mathcal{P}=\left[ \begin{array}{cccc}0&0&0&1\\ 0&0&1&1\\ 0&1&1&1\end{array} \right]. \tag{7}\]  
The first column in the \(\mathcal{P}\) matrix will always be comprised of zeros, since it corresponds to the unfolded (\(U\)(t)) progress curve. As a result, a square and invertible subsection of \(\mathcal{P}\) can be extracted which we denote as **P**. The sub-matrix of Equation 7 is  
\[\textbf{P}=\left[ \begin{array}{cccc}0&0&1\\ 0&1&1\\ 1&1&1\end{array} \right]. \tag{8}\]  
The general form of **P** is a square invertible matrix comprised of ones and zeros, where the columns are linear independent. _A priori_**P** is not known for a given RNA folding reaction since only the experimental progress curves, \(\vec{C_{E}}(t)\), (from a \(\cdot\)OH footprinting experiment) are known. We thus aim to identify **P** and **K** computationally from the \(\cdot\)OH footprinting data. We previously implemented a brute force approach for solving this problem that is described in [9]. In this initial implementation we determined an optimal **K** for all the different possibilities of **P**. For the L-21 _T. thermophila_ group I intron with two intermediates, 28 individual **K** matrices were optimized on separate CPUs corresponding to the 28 different **P** matrices previously enumerated by hand [9]. We used a non-linear least squares optimization that minimizes the difference between the experimental progress curves, \(\vec{C}_{E}(t)\), and the predicted progress curves, \(\vec{C}_{P}(t)\), by adjusting **K** for a given **P** matrix. This process is computationally intensive because we perform non-linear optimization to identify an optimal **K**, requiring multiple solutions of Equation 5 before achieving the convergence criteria. Finally, multiple optimizations are run using different starting values of **K** to correctly identify a global minimum [9]. Together, identifying the winning combination of **P** and **K** can take upwards of 75 CPU days for the _T. thermophila_ group I intron.