Overview.

After the formation of the Cas9:gRNA complex, it must find the target and knock it out. As soon as the complex is formed, it will wander around the nucleus describing a random pathway. During this erratic trajectory, collisions with several regions of DNA will take place. If the union to those regions is thermodynamically balanced and feasible, i.e. regions are complementary enough to the gRNA sequence, the R-loop will take place. This structure results on the hybridization of the Cas9:gRNA complex to a DNA sequence. This implies not only joining the target, but also binding undesired though similar regions, named off-targets.

In this step of the modeling, we used Boltzmann probability distribution and Thermodynamics in order to estimate the probability that the Cas9:gRNA complex finds the target. We also developed an off-target search algorithm based on the transcriptional activity and target-similarity provided by local alignment algorithms.

Complex diffusion.

The process of searching the target among all the genome is named scanning. Thus, we can express the contact rate between Cas9:gRNA complex and any other DNA region:

Where parameters with values indicated in Table 5, are:
# D is the compound diffusivity.
[Cas9:gRNA] is the concentration of the complex.
V is the compartment volume, i.e. the plant nuclear volume.
λ_Cas9 is the characteristic length between the place of production and binding.
Several assumptions were made to consider random three-dimensional diffusion of the complex. Those assumptions were:
# The complex can be considered as a macromolecule with three-dimensional random diffusion around the nucleus.
The net molar flow is presumably equal to zero, as the compartment composition is well-mixed.
As DNA is dispersed among the nucleus, the complex will be almost in permanent contact with it.

Thus, assuming a spherical approach of the Cas9 spatial shape, λ= (V_Cas9)^&frac13, because in the edge of Cas9 there will be DNA ready to hybridize gRNA if possible (2).
Varying the time of measurement and the number of Cas9 and gRNA copies introduced, we can expect different results for this ratio:

Figure 9. Comparison between values obtained for the diffusion ratio of the Cas9:gRNA complex al time t = 4320 minutes (3 days). All curves are relative to the results obtained with 1 copy for Cas9.

Simulations represented in the graphic above let us check that the gain in the rate of contact between the Cas9:gRNA complex, increases approximately in the same way when so does the concentration of the complex. A minimum of 10-15 gene copy number for the gRNA construction is necessary to achieve the plateau of the k_r ratio, independently of the gene copy number for the Cas9 construction. Furthermore, the gain in the number of gene copies encoding Cas9, is the same produced in the random contact ratio.
As it can be inferred from this analysis, achieving enough Cas9:gRNA concentration is critical to stablish contact between the complex and the target. One possibility increase repeatability of the test, minimizing the randomness of the complex diffusion, is to infiltrate the Testing System construction with the gRNA sequence. The expected result was that if the gRNA is transcribed near the target, the aleatory of the three-dimensional diffusion would be minimized. Joining both pieces near to each other, it will be “easier” for the complex to find the target. This suggestion was implemented in wet-lab experiments, showing an increase of the light signal.

Probability of R-loop formation.

In order to knock out our genetic target, it must be hybridized by the gRNA forming the R-loop. This structure provides Cas9 with the necessary stability to cut the DNA strand (2, 15). In order to get the structure, it is necessary that potential targets are complementary enough to the gRNA. Providing gRNA-DNA complementarity means accomplishing the thermodynamic requirements to let the knockout happen.

Estimating the number of targets and off-targets, we can obtain a distribution of the cleavage probability in function of energy needed to cleave each DNA location. Thus, M different energetic states are considered as places where the R-loop could take place, being our Testing System one of that states. This scenario can fit to a Boltzmann distribution (2,5), being the binding probability of the Cas9:gRNA complex with the m-DNA region:
This expression has information about the thermodynamic balance after the R-loop formation. In order to obtain it, we had to obtain previously the free energy increment for each DNA candidate (-ΔG_{complex,target}), and the expected number of those regions (N_target).
Off-target regions were estimated using the off-target search algorithm, getting 1 off-target for Ga20Ox and 5 for TFL. The next section has the explanation of the free energy increment, which ensures the thermodynamic stability of the R-loop.

Off-target search algorithm.

Off-targets are DNA regions where the R-loop could take place because of the high similarity between that region and the target. This means that off targets steal Cas9 and gRNA supposed to knock out on-targets. Most reliable off-target predictions are obtained by experimental results, but in our model we must be able to find off-targets quickly for all possible targets (1,2) so we could not wait for experimental results.
Alternatively, we have developed an off-target search based in transcriptional activities and local alignments between target and off-target candidates. Our proposed strategy was the following one:

The first step was to create a gene library with sequences of the most transcribed genes in Nicotiana benthamiana. There is a clear relation between transcriptional activity and relaxed state of chromatin (13), letting us assume that those genes highly transcribed will be more accessible to the gRNA.
This algorithm is implemented in the Matlab function Nboffsearch.m. The search of potential off-targets for each of our two targets, gave a result of 1 off-target for Ga20Ox and 5 off-targets for the TFL.

Free energy increment ΔGSUBcomplex,targetSB1.

As there is no energy supply catalyzing the R-loop (2), the process described above is a sequence of reactions which in global, must accomplish the thermodynamic law for this kind of processes:

Where the free energy is decomposed in the main stages previously described (2,15):

Each one of the terms from the expression above, is explained in following sections of the Modeling. We determined these parameters for two of the targets contained in our Database: ORYZA SATIVA JAPONICA GROUP GIBBERELLIN 20 OXIDASE 2 (LOC4325003) and CITRUS SINENSIS TERMINAL FLOWER (TFL). Those induce higher grain yield in rice flowering in orange, respectively. Values obtained for each of them were ΔG_{complex,Ga20Ox}= -10,37 kcal/mol and ΔG_complex,TFL =-11,78 kcal/mol.
These values were used to obtain the probability that the Cas9:gRNA complex cleaves the target in the Testing System, performing the desired knockout. To calculate this probability, we had to account for possible off-targets, as they could obtain ΔG_{complex,target} <0, letting the R-loop be formed at those regions as well.
After calculating the ΔG_{complex,offtarget} for off-targets suggested by the algorithm, only the off-target for the Ga20Ox could be considered as a final off-target, as TFL off-targets got free energy increments higher than zero. Therefore, we calculated the P_{complex,Ga20Ox}, obtaining a value of 0,9964.

PAM binding energy..

Once the complex has been formed, it must find the PAM sequence of the target included in the Testing System. The PAM region has between 3 and 5 nucleotides which are recognized by Cas9, enabling the union of the complex to the DNA. Each Cas9 specie binds better to one type of PAM region. In our case, as we are working with the human type Cas9, the PAM sequence will be -NGGN-.
Bearing in mind that breaking one hydrogen bond provides at least an energy supply of 1.2 kcal/mol, the binding of the PAM sequence must give a free energy ΔG_PAM at least of 9 kcal/mol:

Relying on the nucleotides arrangement, ΔG_PAM resulting from the PAM recognition will vary. However, it is possible that Cas9 interacts with a PAM region even though there are mismatches between them (2). In order to know if this affinity for regions with mismatches was significant, we studied the ΔG_PAM obtained for all possible PAM combinations.

White spaces in picture above represent PAM alternatives which do not bind significantly to a -NGG- PAM. The lower binding energy is clearly for PAMs with the structure -NGG-, achieving less than -9kcal/mol. However, there are other alternatives with affinity enough to let the Cas9 bind to them. The potential of these regions as possible off-targets relies also on the Cas9 specie being used, as there are some more “promiscuous” than others. We opted for including this potential off-target PAMs in our off-target search algorithm.
This PAM-energy assignment is implemented in the Matlab function energy_PAM.mat. The input is the target sequence. Comparing the PAM extreme nucleotides of the input sequence with a table containing Cas9-PAM binding energies, it matches the information of the string input with the corresponding energyΔG_PAM . In the particular case of targets implemented in our testing system, the value of this parameter was:
Table 2: Results of ΔG_PAM for targets implemented in Testing System

Cas9:gRNA:DNA hybridization.

Secondly, the release of the energy from PAM binding will be used to hybridize the gRNA and nucleotides from the DNA sequence. In order to know the energy associated to different base pairs, we built a table with the energy increment for all possible matching duplexes. This table provides all information to calculate the ΔΔG_{exchange,gRNA:target} (8,9,10,11), as there should not be mismatches between both of them. In the graphic below it can be appreciated that there are two sources of variability affecting the energetic balance of duplex hybridization. On one hand, there is a clear dependence of the nucleotide, and on the other hand, the type of nucleic acid (RNA or DNA) also affects the free energy increment.

Thus, we had to choose an approach which considered the energy differences between different nucleotides and different nucleic acids. Moreover, the model used to estimate the energetic cost of forming the R-loop, had to consider the distance of base pairs to the PAM region.
This kind of forward move which considers all the mentioned criteria, is known as nearest neighbor model. The meaning of this model is that RNA:DNA union will rely on the context. The energy used to bind a pair of duplexes, uses the energy released by the previous duplex union. In other words, it is assumed that the energy from the k^th hydrogen bond will be used in the reaction of the nearest base pair, k + 1.
Thus, the hybridization of several nucleotides can be represented as a sequence of binding reactions, leading to a global difference between the hydrolysis of DNA:DNA bounds, and the union of gRNA:DNA strands. The term ΔΔG_{exchange gRNA:target} reflects the difference between the free energy used for the pair DNA-DNA hydrolysis and the pair RNA-DNA hybridized.

However, the gRNA may have thermodynamically stable unions with other regions which are not the target. Those DNA regions, named off-targets, may have only few mismatches that slightly affect its affinity towards the gRNA. The half part of the gRNA close to the PAM, is the most determining to form the R-loop. Therefore, if mismatches are placed in the extreme opposite to the PAM, they may will not compromise the off-target knockout.
In a similar way as we did with matching duplexes, our first try was to find more information about energy accounted when mismatches are produced. Nevertheless, there is poor consensus among bibliography, not all possible duplexes have been studied and we neither could determine these energetic values empirically.
In order to solve this, we created a penalty vector which adds a penalty to the match binding energy for each mismatch. The term d_k refers to a weight that decreases as k is increased, with k = 1,2,3…length gRNA (typically 20-23). We have estimated values of those weights using criteria from our Scoring System, as it had been validated comparing to other target searchers available online. Those coefficients are multiplied by the single mismatch average penalty of 0.78 kcal/mol, extracted from bibliography (2). The implementation of the penalties is in the Matlab function weights_exchange.m, and the vector of penalties is represented below

Using this strategy, we studied how would could vary the ΔΔG_{exchange gRNA:target} estimated for a target with different number and positions of mismatches. We implemented the obtention of ΔΔG_{exchange gRNA:target} in the Matlab function energy_exchange.m. Results obtained for our particular targets were:
Table 3: Results of ΔΔG_{exchange gRNA:target} for targets implemented in Testing System.

In order to know more about how do mismatches affect the thermodynamic balance of the gRNA and DNA target hybridization, we calculated values of ΔΔG_{exchange gRNA:target} for rice and orange, varying the position of a single mismatch. Conditions of the simulations are in Table 4.
Table 4: Results of ΔΔG_{exchange gRNA:target} for targets implemented in Testing System.

The results illustrated in the graphic above, show that as expected, the R-loop formation is less likely as it is reduced the distance between a mismatch and the PAM. In general, it seems that our thermodynamics approach emulates well the mechanism of a R-loop formation. With both targets, the average result of changing one nucleotide, is decreased with higher PAM-distance. Mismatches placed downstream the 20 ^th nucleotide, typically result in positive free energy increments, avoiding the RNA:DNA hybridization. This agrees with criteria found in bibliography (2), letting us assume that the penalty system worked well as representation of mismatch effects.
The variability observed in each position is due to the differences between three possible nucleotides. However, there are some atypical results as well which may be caused by unknown sources of variability. For instance, the ΔΔG_{exchange gRNA:target} is overlapped for positions 23 and 10 in rice, while energy difference for the position 23 was supposed to be minor. This could be due to the necessity of training and improving our function, using parameters in the penalty function which are based on empirical evidence.

DNA supercoiling.

Finally, the chromatin state is critical to let the gRNA hybridize the DNA, and some energy can be extra-needed if the DNA is "relaxed", i.e. positively supercoiled. Consequently, regions highly similar may will not be able to join the gRNA because the chromatin could be compressed. Having off-targets means that the binding will be taking place in some unspecified regions. The difference between the initial density of a target σ_I and a non-specific region (σ_NS) will affect to the difference in the free energy needed to untwist the on-target DNA region:

As we could not determine the chromatin state and its evolution during the time that CRISPR/Cas9 was working on the plant, we could not study this parameter.
Nevertheless, information about the chromatin supercoiling has been indirectly introduced in our model. High transcription activities can be synonym of relaxed chromatin (12, 13). Thus, we can assume that supercoiling will not be affecting to our Testing System, since it has the 35S promotor, which is one a constitutive promotor with high activity (4). Moreover, the difference in DNA supercoiling between the target and off-targets, can be considered nearly zero because potential off-targets will be those genes of Nicotiana benthamiana which are highly transcribed.

Parameters.

Table 5: Parameters used calculation of the k_{cleavage TS} parameter.

https://2016.igem.org/Team:Valencia_UPV/Loop

Main remarks

In order to reduce the random influence of three-dimensional diffusion of the complex, we suggested introducing the gRNA and the Testing System constructions one next to the other. Single mismatches positioned downstream the 11^th-10 ^th nucleotide of the gRNA, lead to a positive free energy increment, i.e. they make unable the R-loop between that gRNA and the DNA target.

Conclusions

In this Modeling task, we have developed a mathematical approach of the CRISPR/Cas9 mechanism and its performance in our Testing System. Our model has a set of ODEs representing Cas9:gRNA complex formation, a thermodynamic balance for the R-loop formation, a sequence-alignment algorithm to search off-targets and probabilities estimation to know the reading frame shift after the Cas9 cleavage.
With this in silico model, we have been able to study the behavior of the system under different conditions. For instance, proposals based on the modeling were implemented in wet lab, with some successful results. In particular we determined optimum times of agroinfiltration, and redefined lab protocols for genetic constructions design, as we described below..
Regarding the Cas9:gRNA complex formation, it has been observed in simulations that the fast dynamics of the gRNA causes its fast disappearance. In other words, the gRNA is the limiting reagent. The usual lab protocol considers the simultaneous agroinfiltration of the gRNA and the Cas9. However, in simulations we observed that the production of the complex could be maximized if Cas9 was agroinfiltrated in first place. This is because after 3 days approx. there will be enough Cas9 in the nucleus to sense small quantities of gRNA. This alternative protocol would make the system less susceptible to lower efficiency in the transcription mechanisms of the cell.
The three-dimension diffusion of the complex until it finds the target, describes an aleatory pathway. This implies more variability within the results of the Testing System, making it less robust and reducing its repeatability. In order to reduce this negative impact, we suggested introducing the gRNA and the Testing System constructions one next to the other. Thus, the probability that the complex “randomly” finds the target is increased and there is more control of intrinsic variability sources affecting the result. Wet lab experiments performed under these conditions showed favorable results when the gRNA was agroinfiltrated near to the target.
The analysis of the thermodynamic balance which regulates the R-loop formation, let us estimate that single mismatches positioned downstream the 11^th-10 ^th nucleotide of the gRNA, lead to a positive free energy increment, i.e. they make unable the R-loop between that gRNA and the DNA target. This verifies bibliographic criteria about efficient design of gRNAs.
Another critical point is the chromatin state, which cannot be determined instantaneously nor continuously. We have tried to avoid this lack of information relating the DNA packaging with transcriptional activity. More information about this factor is necessary, letting us further study its relation with the CRISPR/Cas9 efficacy and efficiency.
In conclusion, though lack of data of our testing system under the simulated conditions has been an obstacle for the model validation, the model gives sensible semi-quantitative results. The large amount of time required when working with plants, has limited the quickly generation of data and test of hypothesis extracted from modeling results. Further experimental results in the future will allow to fully validate the model. Nevertheless, results obtained have helped us to determine critical steps of the process, as well as to redesign laboratory protocols.

References

• #1. Bae S, Park J, Kim J. Cas-OFFinder: a fast and versatile algorithm that searches for potential off-target sites of Cas9 RNA-guided endonucleases. Bioinformatics. 2014;30(10):1473-1475.
2. Farasat ISalis H. A Biophysical Model of CRISPR/Cas9 Activity for Rational Design of Genome Editing and Gene Regulation. PLOS Computational Biology. 2016;12(1):e1004724.
3. Chattopadhyay J, Tapaswi P, Mukherjee D. Formation of a regular dissipative structure: a bifurcation and non-linear analysis. Biosystems. 1992;26(4):211-222.
4. Nakabayashi JSasaki A. A mathematical model of the intracellular replication and within host evolution of hepatitis type B virus: Understanding the long time course of chronic hepatitis. Journal of Theoretical Biology. 2011;269(1):318-329.
5. Kalinin MKononogov S. Boltzmann’s Constant, the Energy Meaning of Temperature, and Thermodynamic Irreversibility. Measurement Techniques. 2005;48(7):632-636.
6. Dill K, Ghosh K, Schmit J. Physical limits of cells and proteomes. Proceedings of the National Academy of Sciences. 2011;108(44):17876-17882.
7. Nishimasu H, Ran F, Hsu P, Konermann S, Shehata S, Dohmae N et al. Crystal Structure of Cas9 in Complex with Guide RNA and Target DNA. Cell. 2014;156(5):935-949.
8. SantaLucia, J, Allawi H, Seneviratne P. Improved Nearest-Neighbor Parameters for Predicting DNA Duplex Stability †. Biochemistry. 1996;35(11):3555-3562.
9. Watkins N, Kennelly W, Tsay M, Tuin A, Swenson L, Lee H et al. Thermodynamic contributions of single internal rA⋅dA, rC⋅dC, rG⋅dG and rU⋅ dT mismatches in RNA/DNA duplexes. Nucleic Acids Research. 2010;39(5):1894-1902.
10. Sugimoto N, Nakano M, Nakano S. Thermodynamics−Structure Relationship of Single Mismatches in RNA/DNA Duplexes †. Biochemistry. 2000;39(37):11270-11281.
11. Sugimoto N, Nakano S, Katoh M, Matsumura A, Nakamuta H, Ohmichi T et al. Thermodynamic Parameters To Predict Stability of RNA/DNA Hybrid Duplexes. Biochemistry. 1995;34(35):11211-11216.
12. Deng S, Stein R, Higgins N. Organization of supercoil domains and their reorganization by transcription. Molecular Microbiology. 2005;57(6):1511-1521.
13. Adrian J, Farrona S, Reimer J, Albani M, Coupland G, Turck F. cis-Regulatory Elements and Chromatin State Coordinately Control Temporal and Spatial Expression of FLOWERING LOCUS T in Arabidopsis. THE PLANT CELL ONLINE. 2010;22(5):1425-1440.
14. Ma X, Zhang Q, Zhu Q, Liu W, Chen Y, Qiu R et al. A Robust CRISPR/Cas9 System for Convenient, High-Efficiency Multiplex Genome Editing in Monocot and Dicot Plants. Molecular Plant. 2015;8(8):1274-1284.
15. Rutkauskas M, Sinkunas T, Songailiene I, Tikhomirova M, Siksnys V, Seidel R. Directional R-Loop Formation by the CRISPR-Cas Surveillance Complex Cascade Provides Efficient Off-Target Site Rejection. Cell Reports. 2015;10(9):1534-1543.
16. Moore R, Spinhirne A, Lai M, Preisser S, Li Y, Kang T et al. CRISPR-based self-cleaving mechanism for controllable gene delivery in human cells. Nucleic Acids Research. 2014;43(2):1297-1303.

Downloads