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 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 to undesired regions similar to the target, named off-targets.
In this step of the modeling, we used Boltzmann probability distribution and Thermodynamics in order to estimate the probability that the complex binds to the target, using a search algorithm based on the transcriptional activity and target-similarity provided by local alignment algorithms.
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 λ=VCas91/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.
Varying the time of measurement and the number of Cas9
and gRNA copies introduced, we can expect different
results for this ratio:
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 kr 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
(-ΔGcomplex,target), and the expected number
of those regions (Ntarget).
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
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 ΔGcomplex,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
ΔGcomplex,Ga20Ox= -10,37 kcal/mol and
ΔGcomplex,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 ΔGcomplex,target <0, letting the R-loop be formed at those regions as well.
After calculating the ΔGcomplex,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 Pcomplex,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
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 ΔGPAM at least of 9 kcal/mol.
Relying on the nucleotides arrangement, ΔGPAM 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 ΔGPAM 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
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ΔGPAM .
In the particular case of targets implemented in our
testing system, the value of this parameter is in Table 2.
|Rice - Oryza sativa, Semi-dwarf; higher grain yield.||CGG||-9,600|
|Orange - Citrus sinensis, Induced flowering.||TGG||-9,700|
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 ΔΔGexchange,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
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 kth 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 ΔΔGexchange 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
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 dk 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 ΔΔGexchange gRNA:target estimated for a target with different number and positions of mismatches. We implemented the obtention of ΔΔGexchange gRNA:target in the Matlab function energy_exchange.m. Results obtained for our particular targets were:
|Rice - Oryza sativa, Semi-dwarf; higher grain yield.||-0.7710|
|Orange - Citrus sinensis, Induced flowering.||-2.0785|
In order to know more about how do mismatches affect
the thermodynamic balance of the gRNA and DNA target
hybridization, we calculated values of ΔΔGexchange
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.|
|Orange - Citrus sinensis, Induced flowering|
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.
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 ΔΔGexchange 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
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
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
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.
|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)|
|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|
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 11th-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.