Team:TU Eindhoven/Model
Model overview
dCas gRNA Finder Tool
The dCas9 gRNA Finder tool is used to design very specific gRNAs for the use of dCas9-based detection systems. Searching and filtering of target DNA sequences is done in a target DNA sequence considering a BLAST database of similar coexisting DNA sequences. Two types of detection systems can be used in the program; the paired dCas9-split-NanoLuc system, consisting of two dCas9 subunits, and the single dCas9-BRET system, consisting of only one dCas9 subunit. The tool finds the most specific DNA sequences according to a scoring algorithm and displays them in the GUI. Afterwards, sorting and conversion to gRNA is possible.
Dynamic system
We modeled the dynamics of the interaction between bacteria and phages to determine optimal working conditions for our experiments with phages and bacteria. The model of the dynamic system was created in MATLAB and is based on the article by Santos et. al. “Population Dynamics of a Salmonella Lytic Phage and Its Host: Implications of the Host Bacterial Growth Rate in Modelling” (2014)[1]. We tried to determine model parameters for our proof of concept with E. coli ATCC 11303 and T7 phages were determined by doing several experiments considering the growth of bacteria and phages complemented with bacterial substrate. These parameters had to be estimated by least-squares fitting of the experimental data.
References
- Santos, S. B., Carvalho, C., Azeredo, J., & Ferreira, E. C. (2014). Population dynamics of a Salmonella lytic phage and its host: implications of the host bacterial growth rate in modelling. PloS one, 9(7), e102507.
dCas gRNA Finder Tool
Introduction
We designed two types of dCas9-based detection systems that can bind to DNA fragments with high specificity; the paired dCas9-Split-NanoLuc system, consisting of two dCas9 subunits, and the single dCas9-BRET system, consisting of only one dCas9 subunit. dCas9 uses a guide RNA (gRNA) to find its complementary target DNA sequence. Therefore, for specific detection of a DNA fragment, a highly specific nucleotide sequence must be available in that fragment, which can then be targeted by a designable gRNA. For high specificity to be true, this target nucleotide sequence must also not be available in coexisting DNA fragments in the surrounding environment, called the ‘subject’ sequences. We developed the dCas9 gRNA Finder tool to easily find these very specific sequences in the to be targeted DNA fragment and to design the corresponding gRNA sequences. For naming conveniences, the to be targeted DNA fragment will further on be referred to as the ‘template’ DNA sequence, and the found highly specific nucleotide sequences will be referred to as ‘target’ DNA sequences.
For our proof of concept, we used the lytic T7 phage to replicate in the E. coli ATCC 11303 subtype, and measured DNA concentration before and after addition of bacteria using our sensor system. When designing gRNAs for this process, we inserted the T7 phage genome as template DNA sequence and the E. coli strain B / REL606 genome as subject DNA sequence into our tool. The E. coli ATCC 11303 genome was not available. However, since E. coli strain B / REL606 descends from the same strain, it was assumed that its DNA would be almost identical. When applying our method to e.g. a urinary tract infection, i.a. humane DNA from endothelial cells present in the bladder should also be considered as subject DNA sequences.
The algorithm used by our tool to find target sequences is based on the SSPD algorithm created by the iGEM Peking 2015 team. Depending on the detection system that is used, the general steps for finding target sequences and designing gRNAs are as follows:
- Find all ‘candidate’ target sequences that have a protospacer adjacent motif (PAM) proximal site.
- Select only highly specific candidate target sequences by comparing to the rest of the template sequence and the subject sequence(s).
- Pair left and right candidate target sequences with correct linker length (only for the paired dCas9-split-NanoLuc system).
- Transcript target DNA sequences to gRNA sequences.
Using a BLAST-based two-step filtration, our tool finds the most specific target DNA sequences in the template sequence. The sequences are then sorted on their specificity according to a scoring algorithm and displayed in the tool’s user interface. Sorting on location in the template sequence or on interspace distance can also be performed. Finally, transcription to gRNA can be performed.
A zip file containing the tool and the README file can be found here.
Methods
Find all ‘candidate’ target sequences
dCas9 requires a PAM sequence in the form of 5’-NGG-3’ downstream a target sequence to bind the template DNA, where N is any nucleotide. Python’s regular expression built-in function was used to search for PAM sequences and their corresponding 20 basepair (bp) candidate target sequences. These target sequences are still candidates since they have to be filtered later on. A dCas9 gRNA can bind to the complementary strand of the template strand, or to the template strand itself. Since only the template strand is loaded into the tool, we also search for the reverse complementary PAM sequence 5’-CCN-3’ in the template strand, to accommodate for both cases. A schematic representation of the binding of our paired dCas9-split-NanoLuc system to a template DNA fragment is shown below in Figure 1. Notice that the left guide RNA sequence binding to the complementary strand is identical to the 20 bp target sequence upstream the NGG PAM sequence in the template strand. The right guide RNA binds to the template strand itself, so the 20 bp target sequence in the complementary strand is used for transcription to gRNA. This sequence can be found by complementing the 20 bp sequence downstream the complementary PAM sequence 5’-CCN-3’ in the template strand. To search for NGG sequences, ‘(?<=.{20}).(?=gg)’ was used with the re-function. Searching for CCN sequences was done using ‘(?<=cc).(?=.{20})’ with the re-function.
Figure 1: Schematic representation of the paired dCas9-Split-NanoLuc system binding to a template DNA fragment with gRNAs in PAM-in orientation.
The template DNA is the strand that is loaded into the tool. Note that the 20 bp target sequences are identical to the 20 bp guide RNAs.
Candidate Target Specificity Tests
The specificity of a candidate target sequence is defined as the probability of its corresponding gRNA binding to the target itself and not to other similar non-target sites. This score is quantified by considering both the amount of potential off-target sites and the similarity between the target and its off-targets. Higher scoring targets have a higher specificity, whilst lower scoring targets have a high off-target probability. We adopted a BLAST-based filtration approach (using standalone BLAST) consisting of two steps:
- Filter out candidate targets that have off-targets with an identical 12 bp PAM-proximal sequence.
- Filter out candidate targets with scores lower than the threshold.
1) 12bp PAM-proximal sequence filtration by BLAST
It has been demonstrated in previous research that dCas9 tolerates mismatches better in the PAM-distal region than in the PAM-proximal region. The PAM-proximal bases 8 to 12 largely determine the specificity of the target sequence. Therefore, target sequences that have off-targets with an identical PAM-proximal sequence of 12 bases have a high probability of binding their gRNA at those off-targets. Hence, these sequences are removed from the candidate list. Similar sequences are found by uploading the candidate targets and a sequence to be compared with to BLAST.
When the same target sequence occurs more than once in the template sequence, i.e. there are multiple binding sites, the relative luminescence/BRET will increase much more than if there would only be one binding site. Therefore, to always keep the increase in luminescence/BRET upon detection the same for different template DNA fragments, candidates with overlapping sequences in the self-DNA are also removed. On top of that, in case of the paired dCas9-split-NanoLuc system, these candidates should also be removed considering the later pairing of two candidate sequences. If the gRNAs are also binding at totally different sites, the pair will be formed much less frequently than in a situation where only one target is present. This would require increased concentrations of the sensor system to result in the same amount of paired split-NanoLuc (and bioluminescence level). Therefore, a first filtering round is executed by comparing the candidates to their self-DNA, the template DNA sequence. Then, a second filtering round is executed to remove all candidates with off-targets in the subject DNA.
The remaining candidate target sequences are again uploaded to BLAST to search for remaining potential off-target sites. Sequences found by BLAST that are similar to one of the candidate guide sequences are used in the scoring filter in the second filter.
2) Score filtration
By doing the first filtration step described above, computation in this second step is reduced significantly. For the second filtration step, Zhang Lab score methods were used to give each candidate an off-target effects measurement. For each candidate sequence, the scoring process includes two steps:
- Each of the candidate’s individual off-target sites is given an “individual score”.
- All individual scores are aggregated to an overall score for the given candidate.
The algorithm used to score single off-target sites is shown below:
$$S_{off}=\prod_{p \in M}(1-W(p))\times \frac{1}{\left(\frac{19- \overline{d}}{19}\times4+1\right)}\times \frac{1}{ {n_m}^2}$$
This algorithm consists of three terms;
with
$W=[0;0;0.014;0;0;0.395;0.317;0;0.389;0.079;0.445;0.508;$ $0.613;0.851;0.732;0.828;0.615;0.804;0.685;0.583]$
The above formula takes into account the influence of the mismatches between the off-target and the candidate sequence due to their position (
$$\overline{d}=\frac{\sum_{i=1}^{n_g}M(i+1)-M(i)}{n_g} $$
with the number of gaps ($n_{g}$) defined as:
$$n_g=n_m-1$$
$n_{m}$ represents the number of mismatches in the off-target, which is also considered in
A higher individual score for an off-target indicates a higher similarity to the candidate target, and thus a higher likelihood of gRNA binding to the off-target site. In Figure 2 an example is shown to help readers better understand the calculation of an off-target score. In this case, $M=[3, 7, 12, 19]$ with $W(7)= 0.317$, $W(12)=0.508$, $W(19)=0.583$. Then $\overline{d}=\frac{(7-3)+(12-7)+(19-12)}{3}=5\textstyle\frac{1}{3}$$. This then results in the off-target score:
$$S_{off}=\prod_{p \in [3, 7, 12, 19]}(1-W(p))\times \frac{1}{\left(\frac{19- 5\frac{1}{3}}{19}\times4+1\right)}\times \frac{1}{ {4}^2}=0.0022272111 $$
Figure 2: Schematic illustration of an off-target site. The nucleotides are numbered sequentially from the one most distant to the PAM,
and mismatches in the off-target are emphasized by red dashed lines.
Finally, the candidate sequence score is calculated by aggregating the $S_{off}$ scores:
$$S_{candidate}=\frac{100}{100+\sum_{i=1}^{n_{off}}S_{off}(i)}$$
where $S_{off}(i)$ are the scores for the $n_{off}$ potential off-targets for that candidate sequence. These scores range from 0 to 1 with higher scores indicating higher specificity, i.e. the candidate has less predicted off-target sites and/or off-targets with lower scores. Candidates with individual off-target scores higher than 0.02 or an overall score lower than 0.45 are eliminated due to their lack of specificity. The scores of the remaining candidate target sequences are multiplied by 100% for easier examination.
Pairing left and right candidate sequences with correct interspace distance
For our paired dCas9-split-NanoLuc system, all remaining left and right candidate target sequences are considered for eventual pairing. Two candidates will be paired according to the chosen conformation and appropriate interspace distance. There are four possible conformations; PAM-out, PAM-in, PAM-left and PAM-right (with CCN-NGG, NGG-CCN, CCN-CCN and NGG-NGG PAM pairs respectively). Furthermore, there is the possibility to place the split-NanoLuc proteins at the N- or C-terminus of dCas9. Experimental results show that the combination of C-terminal split-NanoLuc in combination with the PAM-in orientation results in the highest absolute luminescence and in one of the highest fold changes [2]. Since we are working towards a point-of-care device, we decided to prefer high luminescence over a high fold change, and hence, chose this combination of orientation and placement to be the best option for our system. This is also the conformation that is illustrated schematically in Figure 1.
The tool is able to search for all four possible conformations. Since we chose the PAM-in conformation, NGG-CCN PAM pairs will be created. By testing with different interspace distances, the optimal interspace distance turned out to be 70 bp for our paired dCas9-split-NanoLuc system (see Figure 4 on the Results page). The candidates are paired by matching their appropriate positions in the template DNA sequence with the corresponding interspace distance. For the use of the paired system, candidates that cannot be paired at correct distances are eliminated from the list.
Finally, scoring of the pairs is achieved by averaging the scores of the two combined candidates:
$$S_{pair}=\frac{S_{cand1}+S_{cand2}}{2} $$
References
- Hsu, P. D., Scott, D. A., Weinstein, J. A., Ran, F. A., Konermann, S., Agarwala, V., ... & Cradick, T. J. (2013). DNA targeting specificity of RNA-guided Cas9 nucleases. Nature biotechnology, 31(9), 827.
- Zhang, Y., Qian, L., Wei, W., Wang, Y., Wang, B., Lin, P., ... & Cheng, S. (2016). Paired design of dCas9 as a systematic platform for the detection of featured nucleic acid sequences in pathogenic strains. ACS synthetic biology, 6(2), 211-216.
Dynamic system
Introduction
We modeled the dynamics of the interaction between bacteria and phages to determine optimal working conditions for our experiments with phages and bacteria. Phage production needs to be optimized for our system to work efficiently. Since there is a trade-off between workable phage concentrations and speed of the detection, one of the parameters that needs to be determined is the optimal waiting time for phage incubation with the bacteria. Optimal concentrations and dilutions of the bacteria in the to be measured sample also need to be determined when keeping the final application of a point-of-care device in mind. The model of the dynamic system was created in MATLAB and is based on the article by Santos et. al. “Population Dynamics of a Salmonella Lytic Phage and Its Host: Implications of the Host Bacterial Growth Rate in Modelling” (2014)[1]. To determine model parameters for our proof of concept with E. coli ATCC 11303 and T7 phages, we did several experiments considering the growth of bacteria and phages complemented with bacterial substrate. We tried to estimate these parameters by least-squares fitting of the experimental data.
A zip file containing the created MATLAB code can be found here.
MATLAB Model
The dynamic system consists of five delay differential equations considering the bacterial substrate, bacteria and phage concentrations. Phages and bacteria interact with each other in a similar way to predators and prey. Therefore, the mathematical model illustrated here is based on the classical model of Lotka-Volterra. One of the most important features included is temporal heterogeneity, the existence of a time delay in the conversion of bacteria predation in phage progeny. The following five differential equations were used to model the interactions:
\begin{equation} \frac{dS(t)}{dt}=-\alpha\mu(S)X_t(t) \end{equation} \begin{equation} \frac{dX_s(t)}{dt}=\mu(S)X_s(t)-\delta X_s(t)P(t) \end{equation} \begin{equation} \frac{X_i(t)}{dt}=\delta X_s(t)P(t)-\sum_{x=\tau}^{\tau+\rho}B_b(x) \end{equation} \begin{equation} \frac{dX_r(t)}{dt}=\mu(S)X_r(t) \end{equation} \begin{equation} \frac{P(t)}{dt}=\beta\sum_{x=\tau}^{\tau+\rho}B_b(x)-N_i\delta X_s(t)P(t) \end{equation}
Where $S(t)$ is the concentration of bacterial substrate, $X_s(t)$ the concentration of susceptible bacteria, $X_i(t)$ the concentration of infected bacteria, $X_r(t)$ the concentration resistant bacteria and $P(t)$ the concentration of phages at timepoint $t$.
Both susceptible and resistant bacteria grow with a bacterial growth rate $μ(S)$ (Equations 8 and 10). Infected bacteria are assumed not to grow since all their resources are used up for phage production, but they still compete for substrate. Substrate concentration thus decreases due to its usage by the total amount of available bacteria $X_t(t)$, defined as:
\begin{equation} X_t(t)=X_s(t)+X_i(t)=X_r(t) \end{equation}
This decrease is at the rate of the bacterial growth $μ(S)$ multiplied by the amount of substrate needed for a new bacterium ($α$) (Equation 7). Irreversible binding of susceptible bacteria and phages occurs with adsorption rate $δ$, proportional to the product of their concentrations given by the law of mass action. Therefore, both concentrations decrease by this rate while infected bacteria concentration increases with this rate (Equations 8, 9 and 11). Infection of a bacterium can inactivate its remaining receptors [2]. Therefore, we decided to exclude the super-infection phenomenon that may exist when the multiplicity of infection (MOI) is high, and set the amount of phage infections per bacterium ($N_i$) to 1 (Equation 11). Hence, $N_i$ is not considered anymore in further calculations.
The term $\sum_{x=\tau}^{\tau+\rho}B_b(x)$ resembles the delayed release of phages upon bursting of an infected bacterium. Therefore, infected bacterial concentration decreases with this rate and phage concentration increases with this rate, multiplied by the burst size $β$ (Equations 9 and 11). $B_b(x)$ is defined as the number of bacteria that will burst with a latent period of $x$ and is calculated as follows:
\begin{equation} B_b(x)=\delta X_s(t-x)P(t-x)N_d(x) \end{equation}
Here $\delta X_s(t-x)P(t-x)$ is the number of bacteria infected $x$ hours earlier. The factor $N_d(x)$ originates from the fact that not all bacteria have the exact same latent period, and bursts are more likely to be distributed over a certain period. As such, $N_d$ is defined as the proportion of bacteria that will burst with latent period $x$. This proportion follows a normal distribution, where latent period $x$ runs over $τ$ to $τ+ρ$. $τ$ is the latent period at which burst will start to occur and $ρ$ is the rise period, also the width of the distribution. Therefore, this normal distribution has an average of $\overline{x}=τ+ρ/2$ and a standard deviation of $σ=ρ/6$. The calculation of $N_d$ is shown below:
\begin{equation} N_d(x)=\frac{1}{\sigma\sqrt{2\pi}}exp\left(-\frac{1}{2}\left(\frac{x-\bar{x}}{\sigma}\right)^2\right)\Delta x \end{equation}
with step size $∆x= ρ/100$. Finally, while $t\lt x$, $B_b(x)=0$. This makes $\sum_{x=\tau}^{\tau+\rho}B_b(x)$ the sum of 100 burst proportions with each a different latent period $x\in [τ, τ+ρ]$ at timepoint $t$. Due to this factor, the increase of phage concentration is much more smoothened out over a rise period $ρ$ instead of releasing all phages simultaneously at latent period $τ$. Since the bacterial size and phage receptors in the cell wall vary with different growth phases, an equation was also included to enable variation of the adsorption rate:
\begin{equation} \delta\left(\mu(S)\right)=Ae^{\left(\lambda\mu(S)\right)} \end{equation}
where $A$ and $λ$ are parameters determined from experiments. This equation decreases the phage adsorption rate in function of the bacterial growth rate. Considering the equation of $B_b(x)$, the adsorption rate in the delayed term will also be a delayed version with a delayed value of $μ(S):δ(μ(S(t−x)))$.
Experiments and Parameter Fitting
Three experiments were performed to determine the model parameters for our proof of concept with E. coli ATCC 11303 and the T7 phage. All parameters were estimated using MATLAB’s lsqnonlin optimizer. The minimal medium used for culturing was based on M9 medium [3] and consisted of the following: Na2HPO4 6 g/L, KH2PO4 3 g/L, NH4Cl 1 g/L, NaCl 0.5 g/L, MgSO4 0.12 g/L, glucose 5 g/L and 5 mM MgCl2. (Adjusted to pH 7 with NaOH). The full protocols can be found here.
Experiment 1: $μ_{max}$, $K_s$ and $α$
A single-step growth curve of bacteria in minimal medium was created. Pre-inocula of bacteria were cultured in minimal medium containing different bacterial substrate (glucose) concentrations:
- Minimal medium + 0.00 g/L Glucose
- Minimal medium + 0.01 g/L Glucose
- Minimal medium + 0.25 g/L Glucose
- Minimal medium + 0.50 g/L Glucose
- Minimal medium + 1.00 g/L Glucose
- Minimal medium + 2.00 g/L Glucose
- Minimal medium + 4.00 g/L Glucose
- Minimal medium + 8.00 g/L Glucose
Bacterial substrate ($S(t)$ in mg/mL) and bacteria ($X_s(t)$ in CFU/mL) concentrations were measured over time, with time intervals of 20 minutes.
From this experiment, parameters $μ_{max}$, $K_s$ and $α$ were determined. Bacteria multiplication rate $μ(S)$ depends on substrate concentration and is modelled by the Monod kinetics:
\begin{equation} \mu(S)=\frac{\mu_{max}*S(t)}{S(t)+K_m} \end{equation}
From the bacterial concentrations over time, the bacteria multiplication rate over time could be determined:
\begin{equation} \frac{dX_s(t)}{dt}=\mu(s)X_s(t) \longrightarrow X_s(t)=X_{s_0}e^{\left(\mu(S)t\right)} \longrightarrow \mu(S)=ln\left(\frac{X_s(t)}{X_{s_0}}\right)/t \end{equation}
Least squares fitting of the $μ(S)$ and $S(t)$ experimental data to Monod’s model resulted in the parameters $μ_max$ and $K_s$.
The amount of substrate needed for a new bacterium ($α$) was determined by averaging $α$ over the different glucose concentration experiments:
\begin{equation} \alpha=\frac{\Delta S(t)}{\Delta X_s(t)} \end{equation}
Experiment 2: $δ_{exp}$ and $μ_{exp}$
The second experiment followed phage concentrations over a short period time (< latent period) in presence of susceptible bacteria. The experiment is based on the principle that phage adsorption to bacteria decreases the free phage concentration. From these decreasing concentrations, the adsorption rate $δ_{exp}$ had to be determined for a constant bacterial growth $μ_{exp}$. $μ_{exp}$ was determined by measuring the bacterial concentrations before and after the experiment and determining the rate with the same equation (Equation 17) as stated above. Phages and bacteria were co-cultured in fresh minimal medium and phage concentrations were measured at time intervals of one minute by taking samples and adding chloroform to these samples to shut down phage production by bacteria. Afterwards, plaque assays were performed for each sample to determine the decreased phage titer ($P(t)$ in PFU/mL) over time.
$δ_{exp}$ was least squares fitted to the following equation:
\begin{equation} ln\left(\frac{P(t)}{P_0}\right)=-\delta\left(\frac{X_{s_0}}{\mu(S)}\right)\left(e^{\left(\mu(S)t\right)}-1\right) \end{equation}
which was obtained by solving the differential equations below representing phage and susceptible host dynamics respectively:
\begin{equation} \frac{dP(t)}{dt}=-\delta X_s(t)P(t) \end{equation} \begin{equation} \frac{dX_s(t)}{dt}=\mu(S)X_s(t) \end{equation}
Experiment 3: $β$, $τ$, $ρ$ and $λ$
The third experiment follows bacterial substrate, bacteria and phage concentrations over a long timespan to be able to estimate the remaining parameters $β$, $τ$, $ρ$ and $λ$. Again, an overnight pre-inoculum of the bacteria in minimal medium was transferred to fresh liquid minimal medium with 5 g/L glucose. Phage was added immediately afterwards. Measuring phage concentrations over time was done in the same way as described above for experiment 2, but this time with a time interval of 20 minutes. Bacterial substrate and bacteria concentrations were measured parallel at the same time interval.
The remaining parameters had to be estimated by inserting the previously estimated parameters into the full model and least squares fitting the concentration data to the model results.
Results
We attempted to fit the parameters using the experiments described above. Unfortunately, the experiments did not provide realistic and usable results. Hence, we tried fitting all the parameters to the data from experiment 3. This did not turn out to be a great solution either since the fitting process did not succeed. The files used to try the parameter optimization are stored in the ‘Failed’ folder inside the zip file.
Since we were not able to obtain decent parameters that corresponded with literature, we decided to keep the parameters that were given in the model by Santos et. Al. [1]. A plot of the concentrations of bacterial substrate, bacteria and phages, and the bacterial growth rate and adsorption rate was created using the model (Figure 1).
According to literature, the latent period would be around 15 minutes [4, 5, 6]
Figure 1: The concentrations of bacterial substrate, bacteria and phages, and the bacterial growth rate and absorption rate, created using the model.
References
- Santos, S. B., Carvalho, C., Azeredo, J., & Ferreira, E. C. (2014). Population dynamics of a Salmonella lytic phage and its host: implications of the host bacterial growth rate in modelling. PloS one, 9(7), e102507.
- Levin, B. R., Stewart, F. M., & Chao, L. (1977). Resource-limited growth, competition, and predation: a model and experimental studies with bacteria and bacteriophage. The American Naturalist, 111(977), 3-24.
- Sambrook, J., Russell, D. W., Janssen, K., & Argentine, J. (2001). Molecular Cloning: A Laboratory Manual. New York: Cold Spring Harbor Laboratory, Cold Sping Harbor.
- Sillankorva, S., Neubauer, P., & Azeredo, J. (2008). Isolation and characterization of a T7-like lytic phage for Pseudomonas fluorescens. BMC biotechnology, 8(1), 80.
- Storms, Z. J., Teel, M. R., Mercurio, K., & Sauvageau, D. (2019). The Virulence Index: A Metric for Quantitative Analysis of Phage Virulence. bioRxiv, 606350.
- Foster, R. A., & Johnson, F. H. (1951). Influence of urethane and of hydrostatic pressure on the growth of bacteriophages T2, T5, T6, and T7. The Journal of general physiology, 34(5), 529-550.
