A kinetic model predicts SpCas9 activity, improves off-target classification, and reveals the physical basis of targeting fidelity

The S. pyogenes (Sp) Cas9 endonuclease is an important gene-editing tool. SpCas9 is directed to target sites based on complementarity to a complexed single-guide RNA (sgRNA). However, SpCas9-sgRNA also binds and cleaves genomic off-targets with only partial complementarity. To date, we lack the ability to predict cleavage and binding activity quantitatively, and rely on binary classification schemes to identify strong off-targets. We report a quantitative kinetic model that captures the SpCas9-mediated strand-replacement reaction in free-energy terms. The model predicts binding and cleavage activity as a function of time, target, and experimental conditions. Trained and validated on high-throughput bulk-biochemical data, our model predicts the intermediate R-loop state recently observed in single-molecule experiments, as well as the associated conversion rates. Finally, we show that our quantitative activity predictor can be reduced to a binary off-target classifier that outperforms the established state-of-the-art. Our approach is extensible, and can characterize any CRISPR-Cas nuclease – benchmarking natural and future high-fidelity variants against SpCas9; elucidating determinants of CRISPR fidelity; and revealing pathways to increased specificity and efficiency in engineered systems.

Strong off-target sites are identified in silico by a growing set of tools. These tools use bioinformatics 20,21 , machine learning 22,23 , or heuristic 12,14,24,25 approaches to rank genomic sites based on distinctive off-target activity scores. Though such models can identify strong off-targets, they are not quantitative and cannot assess activity on the many lesser off-targets; nor can they predict how activity changes with exposure time and enzyme concentration-even though these parameters are frequently exploited to limit off-target activity in cells 26 .
To implement quantitative activity prediction, Cas9 targeting must be modelled in physical terms. Existing physical models 24,27,28 assume binding equilibration before cleavage, and it remains unclear what predictive power such approaches can ultimately deliver in this non-equilibrium system 29,30 . To account for the nonequilibrium nature of the targeting reaction, we construct a mechanistic model that captures binding and cleavage reactions in kinetic terms. To gain insights into general mechanisms, we train and validate our model on highthroughput datasets that capture both binding and cleavage in bulk experiments 15,31 . Though we restrict our training to offtargets with two or less mismatches, we accurately predict the activities on all more highly mismatched off-targets in the same datasets, as well as those reported in two independent highthroughput datasets 11 .
To reveal the physical basis of Cas9 fidelity on genomic scales, we extract the free-energy landscapes that control PAM binding, strand-replacement, and cleavage on any target. Our characterization of Cas9 supports the notion that observed differences in binding and cleavage activities 32 42 , and our model predicts both its location and its conversion rates.
Though the strengths of our model lies in that it allows us to calculate how (d)Cas9 activity evolves in time under various conditions, we also sought to compare our approach to existing binary off-target classifiers that identify strong off-targets. To this end, we reduce our quantitative activity predictor to a binary offtarget classifier that outperforms the leading tools used today 12,24,28,43 .

Results
The kinetic model. In Fig. 1a we show the reaction pathway that underpins the Cas9 targeting reaction on every target 44 . The reaction starts with Cas9-sgRNA ribonucleoprotein complex exiting the solution state to specifically bind to a 3nt protospacer adjacent motif (PAM) DNA sequence-canonically 5'-NGG-3'via protein-DNA interactions 44,45 . Binding to the PAM sequence (state 0) opens the DNA double helix, and allows the first base of the target sequence to hybridize with the sgRNA 44,45 , forming the first R-loop state (state 1). The DNA double helix further denatures as the RNA-DNA hybrid is extended in the guide-target strand-replacement reaction [46][47][48][49] (state 2-20). The hybrid grows and shrinks in single-nucleotide steps, until it is either reversed and Cas9 dissociates, or it reaches completion at 20 base pairs (bp) in state 20. If the full hybrid is formed, Cas9 can use its HNH and RuvC nuclease domains to cleave both DNA strands 50 .
If we know the conversion rates in Fig. 1a for a particular guide and target, the reaction scheme can be solved to calculate the binding and cleavage probabilities at any time (Methods). Fully parameterizing the model over all guide and target sequences requires the estimation of~10 26 rates. To render parameter estimation tractable, we make four mechanistic-model assumptions: (1) Mismatch positions are more important than mismatch types (e.g. G-G vs. G-A). This can be directly inferred from data 11,15 , and we treat all 12 mismatch types equally. (2) Mismatch energies are determined by local interactions.
The energetic cost of multiple mismatches is taken to be equal to the sum of the energetic costs of the individual mismatches. (3) dCas9 differs from Cas9 only in that dsDNA bond-cleavage catalysis is completely suppressed (k cat = 0); all other rates are taken to be identical 51,52 . (4) All selectivity is governed by the hybrid-bond-reversal rates.
Hybrid-bond-formation rates are treated as equal, independent of complementarity and location.
These assumptions reduce the total number of microscopic parameters to 44 (see Methods): the (concentration dependent) rate of PAM binding from solution (k on ) and the associated freeenergy gain (F 0 ); a single internal forward bond-formation rate (k f ); 20 free energies dictating R-loop progression at the on-target (F 1 ; ; F 20 ); 20 free-energy penalties for mismatches at different R-loop positions (δϵ 1 ; ; δϵ 20 ); and the rate at which the final cleavage reaction is catalyzed for Cas9 (k cat ). Once model parameters are estimated, all possible off-target free energies can be directly calculated using assumptions 1-4 above. In Fig. 1b we illustrate the calculation taking us from the on-target (pink) to the off-target (blue) free-energy landscape with mismatches entering the hybrid at the 3rd and 15th bp. How to translate between free energies and rates is detailed in Methods.
Base-pairing interactions, protein-DNA interactions 52 , and induced conformational changes 50,51,53,54 all contribute to the stability of the Cas9-sgRNA-DNA complex. To account for the varying nature of these interactions, we allow for varying gains and losses in the on-target free-energy landscape as the hybrid is extended. These variable gains and losses allow for the formation of metastable states on the on-target, and constitutes an essential extension of our previous fixed-gain model for RNA-guided nuclease kinetics 30 , as well as of models describing DNA displacement reactions occurring in solution [55][56][57][58] .
Training on binding and cleavage for moderately mismatched targets. We seek to reveal general properties of SpCas9 DNA targeting on genomic scales. To this end, we train and validate our model on data from two highly reproducible bulkbiochemical experiments performed on a large library of moderately to highly mismatched off-targets. The first set 15 (Nucle-aSeq) has 97% correlation between replicated experiments, and estimates the effective cleavage rates (k eff clv ) for a library of offtargets exposed to Cas9-sgRNA for 16 hours. The second set 15,31 (CHAMP) has 94% correlation between replicated experiments, and reports on the effective association constant (K eff A ) over the same library and guide, but this time exposed to dCas9-sgRNA for 10 min. In Methods we detail how the experiments are modeled.
We estimate the model parameters by minimizing the total experimental-error weighted residue between prediction and experiment for off-targets (see Methods) with no more than two mismatches in the NucleaSeq (Fig. 2a-c) and CHAMP ( Fig. 2d-f) experiments. The rates and association constants from different types of mismatches are averaged (see Methods and Supplementary Data 1), and the optimal solution is sought with a Simulated Annealing algorithm 59 (see Methods).
The two training sets differ significantly (Fig. 2, and Supplementary Fig. 1a). Our model still reproduces effective cleavage rates (Fig. 2a-c) and effective association constants ( Fig. 2d-f) with a Pearson correlation of 93% and 98% respectively, and quantitatively captures the difference between binding and cleavage activity. The time and concentration dependence of (d)Cas9 activity can be explored through a dashboard we provide (see Code Availability).
Validation on highly mismatched targets and independent data sets. Apart from the data we use for training (two or less mismatches), the NucleaSeq 15 and CHAMP 15,31 sequence libraries also includes block-mismatched targets with more than two mismatches. In Fig. 3a, b we show that we quantitatively predict effective association constants on these highly mismatched targets at a correlation of 98%. Our method also successfully separates out the single dominating off-target present among highly mismatched targets in the NucleaSeq experiments ( Supplementary  Fig. 1b), resulting in a perfect correlation.
To further validate our model, we test against two data sets from HiTS-FLIP experiments reported in the literature 11 . The first independent validation set records the association rate relative to the on-target, estimated over 1500 seconds of exposure to dCas9-sgRNA at 1 nM concentration (Fig. 3c-e). The second independent validation set records the dissociation rate relative to the on-target, estimated over 1500 seconds following 12 hours of exposure to a saturating dCas9-sgRNA concentration ( Fig. 3f-h). Our model quantitatively captures the relative association rates for all reported targets with 82% correlation (Fig. 3e). For the relative dissociation rates, the correlation is more modest at 46% (Fig. 3h), and the quantitative agreement is lost in some regions (Fig. 3f-h). We still seem to capture the general trends on moderately mismatched targets ( Fig. 3f, g), though our model will never give binding/ dissociation rates above/below that of the on-target, as is reported for some highly mismatched targets (Fig. 3e, h) Physical characterization of SpCas9 and the intermediate R-loop state. As our model parameters carry physical meaning, estimating them from data amounts to characterizing the system in physical terms. For Cas9, it has been experimentally shown that R-loop progression is controlled by an intermediate metastable state on the on-target 42 . We expect this intermediate state to show up as a local minimum in our estimated on-target freeenergy landscape. The free energy of any metastable state will have a strong influence on the observed dynamics, and we expect such energies to be well constrained by the data. We expect barriers between metastable states to be harder to resolve, as the details of barrier regions matter less for the observable dynamics.
We here report 33 near-equivalent optimization runs that all resulted in a residue that fell within 15% of the best solution found (see Supplementary Video 1). In Fig. 4a we plot the resulting on-target free-energy landscapes, with the optimal solution highlighted in pink. As expected, we see metastable states in the on-target free-energy landscape. With Cas9 in solution or PAM-bound, we have a well-defined free-energy minimum where the R-loop is closed (C). The on-target free energy (Fig. 4a) increases substantially when forming the first hybrid bp in state 1, and remains relatively high and poorly constrained up to and including state 8. The energy of state 9-12 are well constrained, and among them we find a second local minimum. We identify these states as belonging to an intermediate NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-28994-2 ARTICLE or reflects the fact that the possible mismatch types vary with position. In Fig. 4c we show the remaining rates needed to predict Cas9 cleavage activity at any target, time, and Cas9-sgRNA concentration (see Methods).

R-loop dynamics captures single-molecule experiments.
The recent direct observation of the R-loop dynamics between metastable states 42 allows us to further test our model against quantitative single-molecule data. To this end, we define a coarse-grained model (Fig. 5a) and calculate the effective rates between metastable states from our microscopic free-energy landscapes (see Methods). In Supplementary Fig. 2 we show that predictions based on our coarse-grained model replicate those of the microscopic model. Using effective rates between metastable states, we can rationalize the broad strokes of Cas9 fidelity by considering a few important examples 42 . For on-targets (Fig. 5b) Mismatches between the target DNA and the sgRNA have differential effects on R-loop propagation depending on position. and future work will have to determine if this is due to a difference in experimental conditions or because our choice of training data is ill-suited to determine the free energies past the intermediate state.
Our model predicts rates on all off-targets, and so extends and refines the long-established rule of thumb that off-target rejection in the PAM proximal seed requires only one mismatch, while offtarget rejection outside the seed region requires multiple mismatches 10 . In particular, our model quantifies the intermediate activity resulting from PAM distal mismatch, and so enables prediction of activity titration.  51 , and considered the time evolution of the occupancy of our metastable R-loop states for two target sequences (Fig. 6). The HNH-domain completes its conformational change within seconds after Cas9-sgRNA binds to on-target DNA 51 , and our model demonstrates a similar behavior for R-loop progression (Fig. 6a) Kinetic modelling improves genome-wide off-target prediction. Current methods 12,14,[20][21][22][23][24][25]28,43 for identifying strong offtargets rank genomic sequences according to various measures of activity. They do not quantitatively predict biochemically measurable parameters, nor do they normally capture changes in conditions or activity over time. Our approach overcomes these limitations, and we do not suggest that these benefits should be abandoned in order to construct a binary off-target classifier. Still, to strengthen the case for including the full non-equilibrium nature of the problem in any Cas9 modelling, we reduce our quantitative kinetic model to a binary classifier (referred to as kinetic classifier) and test how well it performs against three established state-of-the-art off-target predictors: a recent benchmarking of models 28 shows the CRISPRoff classifier to outperform the competition, so we first test against this tool; second, we test against the more recent uCRISPR 24 tool, which is based on hybrid energetics and has not been tested against CRISPRoff; lastly, we test against the Cutting Frequency Determination Fig. 3 Validation on highly mismatched targets and independent HiTS-FLIP data. a Validation data (upper-left triangle) for effective association constant (CHAMP) on block-mismatched targets, and model estimates (lower-right triangle). The two terminal mismatch positions in the block are marked on the axes. b Correlation plot between measured effective association constants and model predictions on block-mismatched targets. c Validation data (triangles) for association rates (HiTS-FLIP data set 11 ) on single-mismatch targets, and model estimates (line). d Validation data (upper-left triangle) for association rates on double-mismatch targets, and model estimates (lower-right triangle). e Correlation plot for all positive association rates, including moderately (1-2 mismatches, dark purple) and highly (3-20 mismatches, light purple) mismatched targets. f Validation data (triangles) for dissociation rates (HiTS-FLIP data set 11 ) on single-mismatch targets, and model estimates (line). The missing mismatch-averaged dissociation rates in the seed are negative. g Validation data (upper-left triangle) for dissociation rates on double-mismatch targets, and model estimates (lower-right triangle). h Correlation plot for all positive dissociation rates, including moderately (1-2 mismatches, dark green) and highly (3-20 mismatches, light green) mismatched targets. Mismatch-averaged rates dominated by negative scores are excluded from the analysis, and all data is averaged over mismatch type (see Methods and Supplementary Data 1). The quoted correlation coefficients are Pearson-correlation coefficients, and correlation plots are displayed with log-scales to show the quantitative agreement also for weak targets. The dashed lines in the correlation plots correspond to perfect quantitative prediction. (CFD) score 12 , since it is a much-used tool for off-target classification. To compare our model against the three selected off-target classifiers, we choose to rank all genomic sites based on cleavage activity in the low enzyme-concentration limit (see Methods). We make our comparison over all canonical PAM sites in the human genome. True positive off-targets are collected from sequencingbased cleavage experiments that used industry-standard sgRNAs and reported multiple off-target cleavage sites [35][36][37][38]40,41,62 (Supplementary Table 1). We tested how well our kinetic model's ranking of activity compares to that of the CFD score 12 , CRISPRoff 28 , and uCRISPR 24 . For each sgRNA, we separately tested the models by using the union (sites found in any experiment) and intersection (sites found in every experiment) sets of the reported off-target sites as true positives. We perform precision-recall (PR) analysis ( Supplementary Fig. 3  b Microscopic free-energy landscape for the on-target exposed to 1 nM (d)Cas9-sgRNA ( Fig. 4a) with coarse-grained states and rates indicated in black. c The calculated (see Methods) coarse-grained forward and backward rates on the on-target. Purple triangles are rates from Ivanov et al. 42 , when available at zero torque. d Microscopic free-energy landscape for an off-target with a mismatch at position 3 (blue), together with the on-target free-energy landscape (pink). Red arrow indicates the free-energy penalty δϵ 3 at the mismatch, and black arrow indicates the resulting shift in barrier height. e The calculated coarse-grained forward and backward rates on an off-target with a mismatch at position 3.
Orange arrow highlights the rate that changed considerably compared to on-target. Purple triangles are rates from Ivanov et al. 42 , when available at zero torque. f Microscopic free-energy landscape for an off-target with a mismatch at position 15 (blue), together with the on-target free-energy landscape (pink). Red arrow indicates free-energy penalty δϵ 15 at the mismatch, and black arrow indicates the resulting shift in barrier height. g The calculated coarsegrained forward and backward rates on an off-target with a mismatch at position 15. Orange arrow highlights the rate that changed considerably compared to on-target. In Fig. 5c, e, and g, central line represents the median, the box plots represent the interquartile range, and whiskers represent the full range among our 33 near equivalent solutions. using receiver-operator characteristics ( Supplementary Fig. 4) since the datasets are highly unbalanced, with many more true negatives than true positives. Figure 7a shows the PR curve when models are tested against the union of all reported off-targets while targeting the HBB gene. As the threshold for what is judged a strong off-target is swept, PR curves display the fraction of predicted off-targets that are found experimentally (precision) against the fraction of experimentally found off-targets that are also predicted (recall). Our kinetic classifier typically produces higher precision for all recalls, outperforming the other classifying schemes for the union set on the HBB gene. More importantly, the kinetic classifier also outperforms the leading off-target predictors for highlymismatched genomic off-targets of other sgRNAs: performing best on the majority of targets in every pairwise matchup on both union (Fig. 7b, c) and intersection (Fig. 7d, e) sets, and irrespectively of if max. F1 or area under the curve (AUC) scores are used.

Discussion
Training our model (Fig. 1) of SpCas9 target activity on moderately mismatched targets, we extract the physical parameters ( Fig. 4) that control activity on any target (Figs. 2 and 3). Going beyond present-day binary off-target classification schemes, we quantitatively predict cleavage and binding activity as a function of both time and SpCas9-sgRNA concentration.
We show that SpCas9's targeting reaction contain an intermediate R-loop state, with both position and conversion rates that agree with single-molecule experiments 42 (Fig. 5). Mismatches affect the dynamics of the R-loop states (Fig. 6) in a manner similarity to how they affect the configurational states of SpCas9's nuclease domains 42,51,53 . Based on this, we lend support to the notion that R-loop formation is tightly coupled to protein conformation-pointing toward the relevant structure-function relation for the most important RNA-guided nuclease in use today.
Though our model captures the abundant low-activity offtargets that are discarded by binary classifiers, we sought to demonstrate the general utility of kinetic modelling by reducing our quantitative activity predictor to a binary classification tool. The resulting kinetic classifier outperforms established state-ofthe-art classification tools on canonical PAM sites in the human genome (Fig. 7).
In a recent study, Jost et al. 5 demonstrated that a series of mismatched guides can be used to titrate gene expression using CRISPRa/CRISPRi. Wildtype SpCas9 can also be (effectively) inactivated with PAM-distal mismatches in the guide 63 . Our model can guide such titration of SpCas9-sgRNA inactivation by careful placement of mismatches. Our approach can also be used to calculate the total off-target activity over a genome, and so inform the design of sgRNAs for novel gene targets.
For simplicity and robustness, we built our model to exclude mismatch type parameters. This allows for extensive training using datasets based on a single guide sequence and off-target DNAs containing up to two mismatches. The limited set of adjustable model parameters (44 in total) and efficient data usage (422 data points used for training) does not seem to limit the model's applicability (Figs. 2, 3, 7). The success of our lowcomplexity model strongly suggest that the path to increased predictive power and therapeutic relevance runs through bottomup modelling of RNA-guided nucleases in kinetic terms.
Taken together, we have shown that our mechanistic and kinetic model gives biophysical insight and quantitative predictive power far beyond the training sets. This predictive power is only expected to increase when including sequence features and allowing for alternative PAM sequences in future modelling efforts. SpCas9 is only one of many RNA-guided nucleases with biotechnological applications, and other CRISPR associated nucleases (such as Cas12a, Cas13 and Cas14) offer a diversified genome-engineering toolkit 15,[64][65][66][67][68][69] . These nucleases can all be characterized with our approach, and it will be especially interesting to compare the free-energy landscape of our SpCas9 benchmark to that of engineered 41,54,70 and natural (e.g. N. meningitides Cas9 71 ) high-fidelity Cas9 variants.

Methods
Modelling of the (d)Cas9 targeting reaction. We consider a single DNA target sequence with a PAM, in contact with (d)Cas9-sgRNA in solution at fixed concentration (Fig. 1a). (d)Cas9-sgRNA binding to the PAM site is assumed to be first order, where [Cas9-sgRNA] is the concentration of active complexes relative to some reference concentration (we use 1 nM). Binding is followed by a Cas9-mediated strand exchange reaction between sgRNA and the DNA. Once a 20 bp hybrid is formed, Cas9 can cleave the DNA, while dCas9 cannot. We model the targeting recognition as a stochastic hopping process along a sequence of states: target unbound (n = −1), PAM bound (n ¼ 0), and strand exchange (n ¼ 1; 2; ; 20). We use the column vector PðtÞ ¼ ðP À1 ðtÞ; ; P 20 ðtÞÞ T to represent the probabilities of being in the various states at time t. The evolution of probabilities is captured by the Master Equation where K is a tri-diagonal rate matrix. Letting k f n be the forward (n ! n þ 1) transition rate, k b n to be the backward (n ! n À 1) transition rate (Fig. 1a), and defining k b À1 ¼ 0,   Fig. 7 Genome-wide off-target classification. a PR curves on the HBB gene using the CFD score (light purple), uCRISPR score (purple), CRISPRoff (dark purple), and our kinetic classifier (green). The precision and recall is calculated over all targets in the genome with a canonical PAM site, taking all experimentally validated off-targets as true positives. b) max. F1 scores for target sites EMX1, FANCF, HBB, RNF2 and VEGFA site 1 using all experimentally identified off-targets as true positives (union set) (Supplementary Fig. 3). c AUC scores for the same target sites and true positives as in Fig. 7b. d max. F1 scores using off-targets identified in all experiments as true positives (intersection set) (Supplementary Fig. 3). e AUC scores for the same target sites and true positives as in Fig. 7d. Matching the models pairwise we can determine which model performs best overall. Using max. F1 scores to count wins on union sets: kinetic:uCRISPR = 4:1; kinetic:CFD = 5:0; kinetic:CRISPRoff = 4:1. Using AUC scores to count wins on union sets: uCRISPR = 5:0; kinetic:CFD = 5:0; kinetic:CRISPRoff = 3:2. Using max. F1 scores to count wins on intersection sets: kinetic:uCRISPR = 2:1; kinetic:CFD = 2:1; kinetic:CRISPRoff = 2:1.
we can give the elements of K as The Master Equation has the formal solution which can be computed numerically, given any set of rates K and initial probabilities Pð0Þ. The above expression, with initial probabilities and rates adjusted to experimental conditions (see below), allows us to capture the full timedependent evolution of the targeting reaction in quantitative terms.
Parameter reduction. Based on the mechanistic-model assumption 1, we average the data over mismatch types (see below), and only keep track of if there is a match or a mismatch at every position. Model assumption 3 means that the model of dCas9 is the same as for Cas9, but with k f 20 ¼ 0. Model assumption 4 implies that To see the implications of model assumption 2, we move to a description in terms of free energies.
Denote the free energy of any state n with F n , and imagine that states n and n À 1 are allowed to mutually equilibrate. Equilibration means that the relative occupancy is described by Boltzmann weights and that there are no net probability currents between the states The above relationships tie rates to free-energy differences through Using n ¼ À1 as the free-energy reference (F À1 ¼ 0 k B T), the assumption that binding is first-order implies is the free energy of the PAM bound state at the reference concentration (1 nM). Mechanistic-model assumption 2 now implies that ΔF 1 ≤ n ≤ 20 only depends on if there is a mismatch at position n or not, and we can write ΔF n ¼ ϵ n ; match ϵ n þ δϵ n mismatch ; n ¼ 1; 20: Here ϵ n is the free-energy increase when extending the hybrid from length n À 1 to length n if the n:th hybrid bp is correctly matched, and δϵ n is the additional energy needed when the bp is incorrectly matched. We can write the backward transition rates as The model is now parameterized it in terms of 41 free energies (F ref 0 , ϵ 1 ; ; ϵ 20 , δϵ 1 ; ; δϵ 20 ) and three forward rates (k ref on , k f , and k cat ).
Predicting NucleaSeq cleavage rates. To produce predications for training and validation, we model experimental setups. To model NucleaSeq data 15 , we use the solution to the Master Equation to calculate the expected cleaved fraction at any complementarity pattern. NucleaSeq is performed by exposing targets to saturating concentrations of Cas9-sgRNA, which we model by setting F 0 ¼ À1000k B T and taking P À1 ð0Þ ¼ 1, P 0 ≤ n ≤ 20 ð0Þ ¼ 0 as initial condition. As done in the original experiment, we record the fraction of DNA that remains uncleaved (∑ 20 n¼À1 P n ðtÞ) at the time points t = 0 s, 12 s, 60 s, 180 s, 600 s, 1800 s, 6000 s, 18000 s, and 60000 s, and fit-out a single effective cleavage rate k eff clv . There is no a priori reason for the uncleaved fraction to follow an exponential decay, but as long as we follow the experimental data-analysis protocol we can use the effective cleavage rates to train and validate our model.
Predicting CHAMP association constants. We model the CHAMP experiments 15,31 by calculating the bound fraction (∑ 20 n¼0 P n ðtÞ) of dCas9-sgRNA after 10 min at concentrations 0.1 nM, 0.3 nM, 1 nM, 3 nM, 10 nM, 30 nM, 100 nM and 300 nM, starting with the probabilities P À1 ð0Þ ¼ 1, P 0 ≤ n ≤ 20 ð0Þ ¼ 0. We use the equilibrium binding fraction to fit out an effective association constant K eff A . Again, there is no a priori reason to believe that this non-equilibrium system will equilibrate within 10 min, but as long as we follow the experimental data-analysis protocol we can use K eff A for training and validation.
Predicting HiTS-FLIP association rates. To predict measured association rates in the HiTS-FLIP experiment 11 , we assume the recorded fluorescence signal to be proportional to our calculated bound fraction of dCas9-sgRNA, when starting with the probabilities P À1 ð0Þ ¼ 1, P 0 ≤ n ≤ 20 ð0Þ ¼ 0. Following the experiments we use linear regression to extract an effective association rate by fitting a straight line to the bound fraction at time points 500 s, 1000 s and 1500 s.
Predicting HiTS-FLIP dissociation rates. To predict measured dissociation rates in the HiTS-FLIP experimen 11 , we again compared the fluorescence signal to our calculated bound fraction of dCas9, starting with the probabilities P À1 ð0Þ ¼ 1, P 0 ≤ n ≤ 20 ð0Þ ¼ 0. We let the protein associate at saturating concentrations for 12 h, and record the resulting occupational probabilities. We then use these probabilities as new initial probabilities, while also letting k on ¼ 0 (½Cas9 À sgRNA ¼ 0) in K, before further evolving the system. This allows us to model complex dissociation in the presence of a high concentration of competitor on-targets in solution. Following the experiments, we fit an exponential decay to our predictions at timepoints 500 s, 1000 s, and 1500 s.
Averaging over mismatch types. Our model does not account for mismatch types, and for training we need to average over all experimentally measured mismatch sequences s consistent with a mismatch pattern p. We expect rates to be proportional to exponentiated transition-state free energies, and association constants to be controlled by exponentiated binding free energies. We therefore choose to perform our mismatch-type averages over the logarithm of rates and association constants, bringing these averages close to averages of energies. For measured quantities m ¼ k eff clv or K ref A , we chose a weighted mismatch-type average hlog 10 m * i p ¼ ∑ s2 sequences with mm pattern p W s log 10 m * s : Here m * s is the measured value for target sequences s. We take the weights to be given by Here δðlog 10 m * s Þ is the experimental error for the logarithm of the measurement at a particular sequence s. This choice of weights minimizes the error-normalized square deviation on the sequence resolved data, if we have complete freedom to set the average for each mismatch pattern. Our model is more constrained then this, but with this weighing our model could-at least in principle-give the best possible approximation of the sequence resolved data. The squared error in the mismatch-type average can be calculated as δ: Cost function. We look to simultaneously optimize our predictions of both effective cleavage rates from NucleaSeq (k eff clv ) and effective dissociation constants from CHAMP (K ref A ). We combine the cost from each experiment by summing log deviations all mm patters used for training w m p ðlog 10 ðm p Þ À hlog 10 m * i p Þ 2 : In the above m p represent the model prediction for the average measured quantity at mismatch pattern p. The weights w m p are chosen so the error-weighted contribution from the on-target, the 20 singly mismatched off-targets, and the 20 Á 19=2 ¼ 190 doubly mismatched off-targets are weighted equally as groups p ¼ on target 1=20; p 2 single mm 1=190; p 2 double mm: Simulated annealing. The Simulated Annealing algorithm 59 is commonly used for high-dimensional optimization problems. We optimize with respect to model parameters F ref 0 , ϵ 1 ; ; ϵ 20 , δϵ 1 ; ; δϵ 20 , log 10 ðk ref on =sÞ, log 10 ðk f =sÞ, and log 10 ðk cat =sÞ. Trial moves are generated by adding a uniform noise of magnitude α to the present value of each model parameter. The process is initiated with a noise strength α ¼ 0:1: In the initiation cycle the temperature is adjusted until we have an acceptance fraction of 40-60% over 1000 trial moves, based on the Metropolis condition. After this initial cycle, the temperatures follow an exponential cooling scheme with a 1% cooling rate (T kþ1 ¼ 0:99T k ). At every temperature, we adjust the noise strength α until an acceptance fraction of 40-60% is reached over 1000 trial moves. Once the desired acceptance fraction is reached, an additional 1000 trial moves are performed to allow the system relax before the next cooling step. Once the temperature has dropped to one percent of its initial value we, apply the stop condition j χ 2 k À χ 2 kÀ1 j ≤ 10 À5 χ 2 kÀ1 : In the above, χ 2 k denotes our cost function averaged over the last 1000 trial moves performed at temperature T k . The results of this optimization is shown in Fig. 4.
Calculating coarse-grained transition rates. First we find the intermediate state on every possible target. As the central-local minimum in free energy (Fig. 4a) can be slightly displaced by mismatches on off-targets, we seek the free-energy minimum n I between R-loop state 7 and 13 for every target. To calculate the effective rates of the coarse-grained model in Fig. 5a, we consider the first passage between metastable states. Take for example the passage from the PAM-bound state (n ¼ 0) to the intermediate state (n ¼ n I ) on a specific target. To calculate the associated first-passage time, we truncate the full system to only include states n ¼ 0; ; n I À 1. We use the rate matrix K PI with elements ðK PI Þ nm ¼ K nm ; 0 ≤ n; m ≤ n I À 1 and k b 0 ¼ 0. With the initial state P PI ð0Þ ¼ ð1; 0; ; 0Þ T we solve the Master Equation, and calculate the first-passage time distribution as Ψ PI ðtÞ ¼ Àð1; ; 1Þ Á ∂ t P PI ðtÞ: The effective transition rate k PI is the inverse of the average first-passage time τ PI , which can be calculated as PI Á P PI ð0Þ: The same process was used to calculate all other rates of directly transitioning between meta-stable states, repeated over every target sequence.
Constructing a binary off-target predictor. We rank all canonical PAM sites in the human genome according to their relative cleavage rate in the low concentration limit. In this limit, the cleavage rate is given by the PAM binding rate times the probability to cleave once the PAM site is bound. As the PAM binding rate is not expected to depend on the sgRNA sequence s, we can rank our offtargets based on the cleavage probability once bound 30

:
Here pðsÞ is the mismatch pattern of sequence s.

Statistics & Reproducibility.
Only experimental data giving physical positive values for mismatch-averaged rates and association constants were included in the correlation analysis. See Supplementary Data 1.
Reporting Summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
The data supporting the findings of this study are available from the corresponding authors upon reasonable request. Mismatch averaged experimental data used for training and validation (Figs. 2 and 3), estimated microscopic parameters (Fig. 4), and genome wide offtarget classification evaluation (Fig. 7b-e), are all provided as Supplementary Data 1.

Code availability
The code enabling quantitative off-target activity prediction for any guide-target pair is available on our GitLab page (https://gitlab.tudelft.nl/depken_group/SpCas9_kinetic _model_dashboard). There you will also find a small dashboard application, allowing time resolved activity predictions given a particular sequence and enzyme concentration. A clone of the repository at publication is also permanently available at https://doi.org/10.5281/ zenodo.5790798. The purpose made optimization code will be made available upon request.