Model - Off-targeting

Introduction to the off-target model

Our project involves employing Cas13a for detecting antibiotic resistance genes. Cas13a acts like a detective that can signal the presence of specific RNA sequences. Its target description is provided by a specific single guide RNA sequence, denoted as the crRNA, that is positioned (loaded) within the protein. Cas13a is activated when the loaded crRNA pairs with a complementary (target) sequence. The crRNA is therefore essential in directing Cas13a towards a predefined target. Besides specific target recognition, proteins from the CRISPR-Cas family may also exhibit off-target activity. While offering possibilities to target related variants of antibiotic resistance genes, this off-target activity could also compromise the reliability of a Cas13a-based detection system. We address this problem by employing a recently proposed kinetic model (Klein et al. 2017) to describe the off-target effects of our detection system for a given crRNA. As such, we are able to rank potential crRNAs - derived from our motif finding algorithm - by their reliability. In addition, these modelling efforts revealed a design flaw in one of the crRNAs, guiding the labteam towards the right crRNA sequences to utilize in our experiments.

The kinetic model we used was designed for proteins similar to Cas9, and required adjustment, as Cas13a has interesting properties that the model was not developed for. Cas13a has a seed region central to the spacer (Abudayyeh et al. 2016), while other Cas protein variants, including Cas9, have a 3’-end seed (Semenova et al. 2011). In addition, for Cas13a, a protospacer flanking site (PFS) – a nucleotide 3’ to the target that should not be guanine - is important, rather than a protospacer adjacent motif (PAM) - a nucleotide sequence next to the non-target strand of for example Cas9 (Shmakov et al. 2015). This can attributed to the fact that Cas13a detects single strand RNA sequences, rather than double stranded DNA sequences. As a result, we set the model’s energy associated to PAM-binding to 0. The process of target binding with a central seed region is illustrated in Figure 1 below.

In the kinetic model, loaded Cas13a can be in a number of states. It can be unbound, have formed $n$ (correct or incorrect) base pairs with a sequence that is being probed, or have cleaved its target substrate, at which point it is also activated. The ratio between the forward and backward rates determines the probability of moving between states and can be determined by the energy difference in the free-energy barrier between going forward and backwards between the state (for each position $n$):

$$\frac{k_{forward}(n)}{k_{backward}(n)} = e^{\Delta n} \tag{1}$$

As mentioned above, functional and structural studies of Cas13a indicate a central seed region (Abudayyeh et al. 2016; Gootenberg et al. 2017; Lui et al. 2017). Therefore, we can schematically represent binding process as depicted in Figure 2). We take the first bound nucleotide to be at position 12 (as seen from the 5’ to the 3’ end of the spacer), in the middle of the central seed region that ranges from nucleotide 9 to 15 (Abudayyeh et al. 2016). We impose that subsequent binding events occur sequentially in both directions.

With these adjustments, we can follow the recipe given in (Klein et al. 2017), we can now calculate the probability to cleave ($P_{clv}$) a target ($s$) with a given guide ($g$) once it is bound:

$$ P_{clv}(s|g) = \frac{1}{1+\sum^n_{n=1} e^{-\Delta T(n)}} \tag{2}$$

With $\Delta T(n)$ the parametrization of the energy landscape: the transition state energy between the first bound state and the n-nucleotides bound state. Lastly, we make the simplifying assumption that the transition to the next state is favourable for correct Watson-Crick base pairing (with an energy gain of $\Delta C$) and unfavourable for mismatches (with an energy penalty of $\Delta I$). The energy landscape is then only dependent on the number of correct ($n_c$) and number of incorrect ($n_i=n-n_c$) base pairs.

$$ \Delta T(n) = -(n-n_c) \Delta I + n_c \Delta C, \qquad n_c\leq n \textrm{ and } n \in [1,N] \tag{3}$$

For convenience, we normalize the $P_{clv}$ through division by the cleavage probability for a matching (on-target) sequence and find R, the cleavage probability relative to the on-target.

$$ \mathcal{R} = \frac{P_{clv} (s|g)}{P_{clv}(g|g)} \tag{4}$$

From a numerical implementation of the model (in Python), we can now determine Cas13a’s tolerance for mismatches between the target and guiding sequence for several model parameters (the energy gain $\Delta C$ and energy penalty $\Delta I$).

Off-target simulations

Although scarce, there are some experimental data illustrating the mismatch tolerance of Cas13a (Abudayyeh et al. 2016, Figure S19 and S21). We observe there is a central seed region between nucleotide 9 and 15. Two consecutive mismatches are only tolerated at the two ends of the spacer, while three consecutive mismatches are only allowed to occur between the nucleotides at the 3’ end of the crRNA spacer and thus the 5’ end of the target. We consider these data, adjust the model parameters and compare the data to the computed tolerance for mismatches using our model. We find a decent parameter fit by eye to be:

$$ \Delta C = 0.25 \tag{5}$$ $$ \Delta I = 2.5 \tag{6}$$

Figure 3 displays the computed mismatch tolerance for single mismatches, two consecutive mismatches and three consecutive mismatches and shows that single mismatches are not tolerated in the central seed region, two consecutive mismatches are tolerated at the ends of the spacer and three consecutive mismatches are only tolerated at the 3’ end of the spacer.

As can be seen from Figure 3, a maximum of approximately 3 mismatches is tolerated by Cas13a. We can use a binomial distribution to compute the probability that a random target and spacer have at most 3 mismatches (thus having 25 or more matches) is low:

$$ P_{\textrm{at least 25 matches}} = \sum ^{28}_{m=25} \frac{28!}{(28-m)!m!} \left(\frac{1}{4}\right)^m \left(\frac{3}{2}\right)^{28-m} = 1.276 \cdot 10^{-12} \tag{7}$$

Having found suitable model parameters for $\Delta C$ and $\Delta I$, we can now employ the model to compute the cleavage probability for different guides and plasmids (or even entire genomes). We started by investigating the relative cleavage probability of our crRNA 1-5 on the tetracycline-resistance conferring (TcR) plasmid that the crRNAs were designed target (Figure 4). We were able to confirm that crRNAs 1-4 target the plasmid at the desired position. However, we also found that crRNA 5 did not meet the design requirements of the PFS: the position 3’ to the target sequence was a G. This finding prevented us from wasting time on testing crRNA 5 in the lab. This is also the reason why we will not consider crRNA 5 in the remainder of the simulations.

By summing the $\mathcal{R}$, the cleavage probability relative to the on-target, over all off-target positions on the plasmid, we compute the $P_{FP}$, the probability of a false positive output, provided that there is a positive output:

$$ P_{FP} = \frac{ \sum _{\textrm{off target}}\mathcal{R} }{\sum_{\textrm{off target}}\mathcal{R} + \sum_{\textrm{on target}}\mathcal{R}} = \frac{1}{1+ \frac{\sum_{on}\mathcal{R}}{\sum_{off}\mathcal{R}}} \tag{8}$$

Computing $P_{FP}$ for crRNA 1-4 (Figure 5) shows that crRNA 3 has the lowest probability of yielding a false-positive result on the TcR plasmid, i.e. being activated by a non-target site on the plasmid: 1 in 100 billion. Interestingly, experiments in the lab show that Cas13a also performs highest collateral cleavage activity when loaded when crRNA 3.

Rewriting the implementation of the model into a command line utility allowed us to use the TU Delft computer clusters for simulating the off-target effects of our crRNAs on the genomes of different relevant organisms (genomes retrieved from NCBI): Escherichia coli BL21 (4.5 mln basepairs (bp)), Staphylococcus aureus (2.8 mln bp), Staphylococcus epidermis (2.5 mln bp) and chromosome 27 of the domestic cow (Bos taurus, 45 mln bp). E. coli BL21 is a commonly used lab strain, which was the strain containing the tetracyclin resistance-conferring (TcR) plasmid that we designed the crRNA for. S. aureus and S. epidermis are common pathogens in mastitis, while genetic material from the cow itself may also be present when detecting sequences in cow samples, for example milk. Scanning the entire genome allows to investigate if tiny local off-target effects together influence the reliability of the output signal.

In Figure 6, we present $P_{FP}$ for our crRNA sequences for the different genomes/chromosome. In general, we find off-target cleavage probabilities below 1 in 100.000 for the bacterial genomes. As expected, the $P_{FP}$ is generally larger for the Bos taurus chromosome 27, which is approximately tenfold larger than the bacterial genomes, but reaches a maximum of 1 in 1000. As such, we conclude that the probability of a false output is low and likely does not compromise the usability of the Cas13a-mediated detection system.

The simulations so far indicate the detection system to be reliable, i.e. have a low probability of yielding a false-positive output. However, insufficient data are available to accurately determine the model fit parameters. In addition, experimental validation is lacking. As such, it is possible that experimental work proves Cas13a to have increased off-target activity. We investigated what the effect of this would be on the reliability of the system, and subsequently propose methods to experimentally reduce the off-target probability.

We repeated the simulation for a higher energy gain for matches ($\Delta C=0.30$) and lower energy penalty for mismatches ($\Delta I=0.80$). In this situation, we find the off-target probability to be much higher than in our original simulations: $P_{FP}$ is 10 to 800 (Figure 7). It is now more intuitive to interpret $P_{FP}$ as the total number of false positives per true positive. Given a positive output, the probability of it being false is higher than 90% in all simulations. This situation would greatly compromise the reliability of the detection system, so we set out to investigate possible solutions.

In our simulations, we assume equal representation of all genomic sequences. In reality, only part of the sequences are transcribed and available in RNA, and of course in different quantities. We can correct for these effects by applying a mask over the genome which takes into account the amount of RNA transcripts for each position. Moreover, our simple sample preparation method for gene amplification (RPA) and in vitro transcription are able to specifically amplify a region of interest at least 452 times (notebook of sample prep 15-08-2017). This amplification changes the ratio between true and false positive output probabilities (by the same factor of 452) in favour of on-targeting (Figure 8). We can incorporate this amplification factor in a generalized version of $P_{FP}$: $P_{FP}'$ and we call the amplification factor α:

$$ P_{FP}' = \frac{1}{1+\frac{\alpha \sum_{on}\mathcal{R}}{\sum_{off}\mathcal{R}}} \tag{9}$$

From our experiments, we estimate α$=452$. In the case of no amplification, we can set α$=1$. The amplification step improves the reliability of the system from 10% up to 40-95%, even in the case of higher mismatch tolerance ($\Delta C=0.30$,$\Delta I=0.85$)

As such, gene amplification and sensitivity not only increase the sensitivity of the detection method, these two steps in sample preparation also reduce the probability of false positive results and thereby increase the reliability of the detection system. We can now simulate $P_{FP}^\textrm{'}$ for different ratios between on- and off-target cleavage probabilities, and determine the amplification factor required to achieve reliable detection above a desired threshold $\beta$ (Figure 9). Analytically, this amplification factor can be calculated as follows by setting $P_{FP threshold}'=\frac{1}{\beta}$:

$$\alpha_{required} = \frac{\beta}{1-\beta} \frac{\sum_{off}\mathcal{R}}{\sum_{on}\mathcal{R}} \tag{10}$$

We can see the reliability of the system increase as a function of α for different rations $\left( \frac{\sum_{on}\mathcal{R}}{\sum_{off}\mathcal{R}}\right)$ and can read-off the amplification factor required to achieve an input that is 95% reliable (i.e. the probability of a true positive given a positive output).

Collaboration iGEM LMU-TUM Munich

Both the iGEM TU Delft and iGEM LMU-TUM Munich teams realized that optimal crRNA is a prerequisite in utilizing Cas13a for their respective detection purposes. In addition to sensitively detecting a sequence of interest, modeling the off-target effects of a chosen crRNA is important in assessing the reliability of the system, as it provides information about the probability of a false positive output. Both teams took a different approach: we employed a kinetic model (Klein et al. 2017) while iGEM LMU-TUM Munich applied an alignment algorithm.

Biophysical modeling has the advantage of providing a unified targeting framework that is able to take into account tiny (local) off-target effects that together may add to affect the reliability of the system. However, it is computationally intensive as compared to homology search and requires the development of the software, which is already available from NCBI for homology search.

We ran 4 off-target simulations for 2 different crRNAs spacer sequences (crRNA1-3 and crRNA4) designed by iGEM LMU-TUM Munich for targeting the 16s ribosomal RNA of two different bacterial species (E. coli DH5α and B. subtilis). The on and off-target activities are modeled for sequences on all possible frames on the genome (Figure 10). From the simulation results, we deduce that the crRNAs specifically target the different bacterial species.

We establish that the two designed crRNA sequences target their designated bacterial species specifically. crRNA 1-3 contains 7 on-target hits on the E. coli DH5α genome, while crRNA 4 contains 10 on-target hits on the B. subtilis genome. From this it can be concluded conclude that the designed crRNA sequences allow differentiation between the (Gram-negative) E. Coli and (Gram-positive) B. subtilis species. Arguably, designing one crRNA that targets both strains could be beneficial in distinguishing between for example a viral and bacterial infection by requiring less experiments, however a Gram-positive or negative distinction could is already very helpful.

We also compute $P_{FP}$ for the four simulations (Figure 11). For all simulations, the probability of acquiring a false-positive result was found to be low: in the order of one in a million times. Moreover, having multiple on-target sites on the genome of the bacterial species to be detected decreases $P_{FP}$ on the targeted genome. The $P_{FP}$ on the non-targeted genome, however is unaffected. This implies that having more on-target sites intuitively also increases the sensitivity of the assay. However, the off-target probability on non-target species remains unaffected by this.

We conclude that the design of the crRNAs is reliable as only minute off-target effects in the order of one in a million are observed. In addition, we gain a deeper insight into the effect of having multiple on-targets in the genome: this increases the sensitivity of the detection method, but does not affect the reliability. As such, our tool assesses the design of the crRNAs by iGEM LMU-TUM Munich to be clever and to allow reliable detection. We thank iGEM LMU-TUM Munich for the fun collaboration and look forward to meeting them again at the Giant Jamboree.

Conclusions

We employed a recently proposed kinetic model to assess the reliability of our crRNAs in the Cas13a-mediated detection system. In general, we find the predicted false positive probability (given a positive output) to be low (as low as 1 in billion). Moreover, the simulation unveiled a design flaw in crRNA 5 that prevented us from wasting time on testing it in the lab. Subsequently, we simulated the off-target activity of all our crRNAs on the genomes of relevant organisms: Escherichia coli BL21, Staphylococcus aureus, Staphylococcus epidermis and chromosome 27 of the domestic cow (Bos taurus). As expected, the false positive probability increases for the presence of larger genome sequences being present. We repeated the simulations for parameter that correspond to a higher mismatch tolerance and observed that this resulted in false positive probabilities exceeding 90%. This situation would greatly compromise the reliability of the detection system, so we set out to investigate possible solutions. In order to do so, we re-evaluate the model’s assumptions and find that if we incorporate our sample preparation protocol, amplifying and in vitro transcribing the gene, the probability of false positives is reduced accordingly. In addition, we propose a method that allows calculation of the required amplification factor α for achieving a desired reliability in the output.