Model and assumptions

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.

The picture below represents a spherical approach of the Cas9 shape, being λ=V_Cas9^1/3. We can consider that next to the edge of Cas9 there will be DNA ready to hybridize gRNA if possible (2). Assumptions made to consider this situation are:

• 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.

Simulations

Varying the time of measurement and the number of Cas9 and gRNA copies introduced, we can expect different results for this ratio:

Figure 10. 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 to 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 ΔG_{complex,target}.

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.

Figure 11. Energies obtained by combining all possible RNA and DNA pairs with different PAM sequences. Units are kcal/mol.

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 is in Table 2.

Target	PAM	ΔGPAM
Rice - Oryza sativa, Semi-dwarf; higher grain yield.	CGG	-9,600
Orange - Citrus sinensis, Induced flowering.	TGG	-9,700

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 necessary 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.

Figure 12. Comparison between the energy exchange resulting from hybridizations with all possible nucleotides of RNA or DNA

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 in the equation above, 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, collected in a vector which is represented on the right, 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.

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:

Target	ΔΔGexchange gRNA:target
Rice - Oryza sativa, Semi-dwarf; higher grain yield.	-0.7710
Orange - Citrus sinensis, Induced flowering.	-2.0785

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.

Original nucleotide	New nucleotide	position	ΔΔGexchange gRNA:target
Rice - Oryza sativa, Semi-dwarf; higher grain yield.
G	A	23	-0,563
G	T	23	-0,603
G	C	23	-0,703
G	A	10	0,369
G	T	10	0,169
G	C	10	-0,381
C	A	7	4,989
C	G	7	2,389
C	T	7	4,589
Orange - Citrus sinensis, Induced flowering
G	A	23	-1,8805
G	T	23	-1,8905
G	C	23	-1,9705
T	A	12	-0,8885
T	G	12	-0,3885
T	C	12	-1,3385
C	A	4	2,6815
C	G	4	-0,3185
C	T	4	1,6815

Table 4: Results of ΔΔG_{exchange gRNA:target} for targets implemented in Testing System.

Results illustrated in the graphic below, 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.

Figure 13. Estimation of total energy exchange produced between the gRNA and the DNA target, using examples of TFL and Ga20ox.

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.

Since high transcription activities can be synonym of relaxed chromatin (12, 13), we can assume that supercoiling will not be affecting to our Testing System, as it has the 35S promotor (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

Parameter	Value	Source
D	2700 µm2/min	Reference (6)
λ	0,015 µm	Reference (7)
V	14140 µm3	Waterloo iGEM team 2015
[Cas9:gRNA](t)	t = 3 days	Model estimated.
ΔΔG single mismatch penalty	0.078 kcal/mol	Reference (2)
kr	0.0172⋅[Cas9:gRNA]	Model estimated
kc	0,48 min-1	Reference (2)
kunbind	300 min-1	Reference (2)
kB	0.0019872041 kcal/(mol⋅K)	Reference (5)
N(OFF-TARGETS)	1 for Ga20Ox	Model estimated
T	297 K	Experiment conditions
P(complex,Ga20Ox)	0,9964	Model estimated

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

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.