Model

Bioinformatic Analysis

The first step to lay the foundations of our project was based on the searching and validation of miRNA biomarkers highly expressed in patients' blood samples compared to control samples. To achieve this goal, we applied a bioinformatic analysis on the dataset with “GSE137140” accessed by Gene Expression Omnibus database. Utilizing R Studio and miRNA Differential Expression Analysis (miRNA DEA) under the R package limma (Ritchie et al., 2015), 1931 miRNAs were significantly over-expressed in lung cancer samples compared to controls fulfilling the following criteria: a) \( P_{\text{val}} \lt 0.05 \), b) \( \log(\text{FC}) \gt 2 \) (log Fold Change). Furthermore, based on the model introduced by Keisuke Asakura and colleagues the same dataset revealed that miR-17-3p is the best single miRNA biomarker for the detection of lung cancer in blood samples with 93.3% sensitivity and 88.5% specificity (Asakura et al., 2020). miR-17-3p displayed a high score of differential expression in our analysis as well, showing \( \log(\text{FC}) = 6.397699 \) (which is almost 84 folds overexpressed in lung cancer samples compared to controls), and \( P_{\text{val}} = 0 \), a strong statistical significance. Therefore, we chose to study our diagnostic system with miR-17-3p and miR-17-5p as well, since previous reports revealed that both miR-17 strands work synergically (Borzi et al., 2021) and mentioned the considerable value of miR-17-5p biomarker in lung cancer detection in early stages (Journal of Cell Science, 2013). Further details about the bioinformatic analysis are displayed in the software section.

Design and evaluation of CRISPR - Cas13a system

This part is crucial for the proper function of our system. We aim to design sequences that will exhibit good thermodynamic behavior and will function in an appropriate manner in order to bind to LbuCas13a protein and also perform detection of our candidate biomarkers.

Design and Analysis of crRNA sequences

In Figure 1 is demonstrated the 2 main parts of every guide RNA or crRNA sequence. The first part is called Repeat part and consists of the nucleotides that will bind to LbuCas13a protein, and it is a specific sequence. The second part is called Detection part and is responsible of detecting the target RNA that is designed to interact with. The length of the detect sequence contains 20 nucleotides that is complementary to the target RNA sequence (Wang et al., 2020).

For the purposes of our research project the biomarkers that were chosen were miRNA-17-5p and miRNA-17-3p. The sequences of these miRNA are demonstrated in Table 1. The information about the sequences were obtained from miR Base.

Based upon these sequences we can design crRNA candidates for detecting these miRNAs. The sequence of every candidate consists of, as described before, 2 parts. The repeat part is the same for every sequence and the detection part consist of 20 nucleotides that are complementary with these biomarkers.

After the design of the candidate sequences, we need to discern the best candidates to select for the experimental part. The criteria for the evaluation of these candidates were: 1) Minimum Free Energy, 2) GC content, 3) Linearity of the detection part, 4) Binding Energy and 5) Perfect Matches (Chau et al., 2020).

The first 3 criteria are based on the secondary structure of the crRNA candidates.

• Minimum Free Energy (MEF) is demonstrating the structures' stability of the sequence. The more negative the MEF is, the more stable will the structure be. The MEF is calculated with software that evaluates the best folding conformation of the sequence and then rank every conformation with a given energy. The energy of every structure is described with the term Delta G (Gibbs free energy). The software that was used for the calculation of the best folding conformation for its sequence was the RNAfold from ViennaRNA (Lorenz et al., 2011).
• GC content is ratio of the Guanine and Cytokine that are present in the sequence. The higher the ratio, the more stable the sequence will be. This phenomenon is based on the different chemical groups that are present in every nucleotide resulting stronger hydrogen bonds in comparison with other nucleotides.
• Linearity Of Detection Part (LODP) shows the secondary structure of the detection part of the crRNA candidates. The secondary structure is evaluated through the energy term of Gibbs free energy of this part of the sequence. Higher the energy more linear will this part be. The linear structure of the detection part of crRNA is highly associated with the proper function of the CRISPR/Cas13a system.

The final 2 criteria are based on the interaction between the crRNA sequences and the according target miRNA (Chau et al., 2020).

• Binding energy is the energy term of the difference between the Gibbs free energy of the complex RNA and the free energy of crRNA and miRNA target. The binding energy is calculated using the secondary structure of the sequences. The software used is the RNAcofold from ViennaRNA.
• Perfect Matches is a measure of evaluation of the RNA complex that is formed by crRNA and miRNA target sequences. The secondary structure of the complex is calculated through RNAduplex, and perfect matches is the ratio of the nucleotides that are bind in detection part of the crRNA sequences.

In Table 4 and 5 are demonstrated the results of the criteria by the ViennaRNA software.

By the demonstrated results from Table 4, it can be concluded that for the crRNA candidates from miR 17-3p, all the candidate sequences have the same GC content, LODP and Perfect Matches. So, the difference for these sequences is about the MEF and the Binding Energy. It can be discerned that the difference between the values of MEF are smaller in comparison with the difference between the values of Binding Energies. Consequently, the ideal candidate sequence for the miR 17-3p is the cr_17_3p_2. Based on the same criteria the ideal candidate for the miR 17-5p is cr_17_5p_1.

After selecting the candidate sequence for the target miRNAs, the next step was to design 2 new generation of sequences. The second generation of the sequences was with the introduction of an artificial mismatch in the detection part of the sequences. The artificial mismatch in the detection sequences was introduced in the position 14 in the detection part for every sequence. The replacement nucleotide for this position was the same nucleotide with the target miRNA sequence, to also prevent the wobble GU base pairs. The third generation of sequences involve the introduction of extra nucleotides in the end of the detection parts. The extra nucleotides were introduced to achieve an extra hairpin loop and to minimize the off-target effects. In the Table 6 and 7 are shown the sequences with the different designs (Ke et al., 2021).

After the design of these different generations of designs, we proceed for the analysis of the thermodynamic behavior of these sequences. In Table 8 and 9 are demonstrated the results of the analysis for the targets miR 17-3p and miR 17-5p respectively. The analysis was based on the 5 criteria as the initial one.

In Figures 2 and 3 are demonstrated the figures of the folding sequence. In the left column are shown the folding sequences alone that were implemented with the software RNAfold. In the right column are shown the folding sequence with the target miRNA, green for the crRNA sequences and red for miRNA targets. These figures were obtained through the RNAcofold software (Lorenz et al., 2011).

Standard

Mismatch

Hairpin Loop

Standard

Mismatch

Hairpin Loop

As can be seen from the results, the extra hairpin loop at the end of the detection part of the 3rd generation crRNA sequences (cr_hs_17) are resulting in lower MEF. This reduction of the energy comes from the extra base pairs making the structure more stable in comparison to other 2 generations. From the results of the Binding Energy evaluations can be seen that the 3rd generation exhibits the lowest binding energy. This is expected since the binding energy is the difference between the free energy of the complex and the free energy of the sequences alone. Since the cr_hs have lower energy and since they don't form more hydrogen bonds, the binding energy is lower in comparison with the other sequences. The other notable characteristic of the results are the perfect matches. The ratio of the designed perfect matches changes only for the artificial mismatch sequences, since the mismatch is making the formation of the hydrogen bond between the nucleotides impossible. Finally, the Linearity of Detection Part is decreasing for the second and the third generation. This can be explained through the different structures that are adopted from the different generations of sequences, since the structure with the minimum energy is preferred.

The final step of this analysis was to determine if it was possible to predict the Level of Detection (LOD) for every sequence. The modelling approach that was selected was through the concentration dependency plots. Concentration dependency plot is a plot that depicts the concentration of the complex based on the concentration of the reactants RNA elements at equilibrium state. In Figure 4 are demonstrated the plots for every crRNA sequence with the according target miRNA.

For every diagram the green line depicts the relative concentration of the complex and all the other lines are covered from the yellow. By seeing these diagrams, we can conclude that in equilibrium all the crRNA sequences will form complexes with the according miRNA targets. The number of the complex will be at their theoretical maximum values since the relative concentration can exceed the value of 0.5. It can be easily understood that utilizing the ViennaRNA software it is not possible to predict the LOD. This is due to the mathematical model that is used for the calculation of the concentration for every RNA element. The model is based on the free energy for every possible state (AA, BB, AB, A, B with A to be the crRNA sequence and B the miRNA). The free energies of the complexes are far greater than the energies of the other 4 states for every crRNA sequence. Consequently, the preferable state of the system will be the formation of the complex in every concentration.

crRNA/miRNA/Protein interactions - Molecular Docking

The next step of the analysis of the system is the evaluation of the interactions between the crRNA sequences and the formed complexes with the protein LbuCas13a. The process of evaluation can be studied through molecular docking (Wang et al., 2020).

Molecular docking is a method which aims to predict the preferable orientation of one molecule to a second when one molecule is chosen to be the receptor and the other the ligand. The purpose of this procedure is to predict the formation of a stable complex. Molecular docking evaluates all the possible conformation of the receptor and the ligand with a scoring function that is associated with energy terms. The purpose of molecular docking software in general is to predict the recognition and interaction between the molecules. The types of docking are based on shape complementary methods or ab initio methods to conclude which conformation is the most stable. For our purpose 2 online web- tools were used, the HDOCK server and the ZDOCK server. The HDOCK server is based on a hybrid docking algorithm of template based-docking and free docking and on the other hand the ZDOCK server is based on a rigid docking algorithm (Zhang Di Yumeng Yan, 2017, Kasprzak et al., 2020).

Before the presentation of the docking results and methodology is appropriate to describe the conformation of the LbuCas13a - crRNA system. The protein LbuCas13a consists of 2 main lobes, the Recognition Lobe (REC) and the Nuclease Lobe (NUC). The Recognition Lobe contains 2 domains, the NTD and Helical-1 domain. The NUC Lobe contains 4 domains, HEPN1-I, Helical-2, HEPN1-II, Linker and HEPN2. From the X-ray data of complex LbuCas13a/crRNA can be seen that the crRNA sequence is bind to the REC Lobe and when the target RNA is present near the complex, the complex changes its conformation and the complex between the crRNA and the target RNA is situated in the NUC Lobe. After this change of the conformation, the protein is activated and the trans- catalytical activity begins (Liu et al., 2017).

The purpose of the docking analysis is to evaluate if the binding conformation that is predicted with the docking algorithms is consistent with the proposed binding conformation of the literature. The main steps of the docking analysis process are:

1. PDB structure of the protein
2. PDB structure of the RNA sequences
3. Fix of the protein missing residues and atoms
4. Docking algorithm
5. Docking Results evaluation.

For the purposes of acquiring the pdb file of the protein, it can be achieved through 2 different ways. The first from Protein Data Bank. The Protein Data Bank ID of the 2 states of the protein (with only the crRNA and with the complex crRNA and target RNA) are 5XWY and 5XWP respectively. The second one is for the AlphaFold web server, where the residues of the protein chain can be used as an input and the structure can be obtained with the machine-learning algorithm of the AlphaFold server. To achieve the exact same structure as the literature, the first one was used.

For the acquisition of the PDB files of RNA sequences the RNAComposer webserver was used. As inputs utilized the nucleotides of the sequences and the secondary structure of the MEF for the energy minimization purposes with dot-bracket notation.

In the following figures are shown the protein structure with its lobes with different colors (Blue - REC Lobe, Yellow - NUC Lobe). In the next figure are shown the different generation of crRNA sequences for miR 17-3p, with red color is the Repeat Part of every sequence and with green is the Detection Part.

Standard

Mismatch

Hairpin Loop

After the acquisition of the pdb file of the molecules (crRNA sequence and protein) we need to prepare the receptor molecule for the docking procedure. The preparation of receptor molecules - protein of interest - contains the repair of the missing residues and the missing atoms of the molecule.

After the preparation of the protein structure the docking webtools can be used. We utilized 2 different webtools for docking purposes: the HDOCK and ZDOCK server. These two tools have different approaches for predicting the best pose of the receptor and ligand complex. ZDOCK utilizes a rigid docking algorithm and the HDOCK server utilizes a hybrid algorithm of template based docking and free docking.

In the following figure are demonstrated the results of the energies of the complexes of the different sequences with the protein that was calculated from the HDOCK and ZDOCK server. The best pose for every sequence was selected to evaluate the results of these 2 different methods of docking.

The first thing that needs to be noted is that the calculated energies differ. The difference between these docking algorithms is on the scoring function that is utilized by them. Finally, it can be shown that the shape complementary rigid method that is used for ZDOCK can lead to unwanted poses such as in Figure 7. This happens because of the shape of the pdb file and because they don't match exactly. For that reason, we decided to continue with the HDOCK server.

Since, in some cases is best the usage of more than 1 pose of the complex, 2 states were evaluated for every interaction on the HDOCK server. The first was the evaluation of the best model of the HDOCK. The second one was the evaluation of the combination of the best three models. After selecting the method of studying the docking results, we proceed to evaluate the interactions between the sequences and the protein and mainly the hydrogen bonds formed between the molecules.

In the Figure 8 below can be seen the results of the Docking score and Ligand RMSD for the best three models for every sequence. It can be understood that the combination of more than one model are unnecessary since the information about the Docking Score and Ligand do not change. In Figures 9 and 10 are demonstrated the intermolecular interactions in the 2 different Lobes of the Protein and the 2 different parts of every sequence. Also are demonstrated the hydrogen bonds that are formed in the 2 Lobes with the BIOVIA software.

Standard

Mismatch

Hairpin Loop

Standard

Mismatch

Hairpin Loop

Standard

Mismatch

Hairpin Loop

Standard

Mismatch

Hairpin Loop

Standard

Mismatch

Hairpin Loop

Standard

Mismatch

Hairpin Loop

It can be shown that most of the interactions between the sequences and the protein take place on the NUC Lobe. This is not in agreement with the literature, but this can be caused by the different shape of the RNA chain in comparison with the RNA chain of the X-ray data. Also, most of the interactions are between Repeat part of crRNA sequences in comparison with the detect part which agrees with the literature. In general, it can be concluded that the usage of the 3 best models, do not improve the prediction of the actual binding site of the crRNA sequences. In Figure 11 are demonstrated the ratio of the off-target interactions for every sequence with the protein. This evaluation was based on the interactions between the PDB protein LbuCas13a 5XWY with the crRNA from crystallography data. The mathematical formula used for the calculation of the ratio was:

It can be concluded that values of this ratio near to 0, then the better the docking model will be. From the Figure the docking models failed to interpret the interactions of the sequences with the REC Lobe and all the interactions made were with the NUC Lobe.

The final step of the docking analysis is the tertiary system docking. In the procedure we need to obtain the pdb file of the complexes crRNAs and target miRNA and the pdb file of the LbuCas13a with the different conformation. The pdb file of the protein is obtained through Protein Data Bank with ID 5XWP. The pdb file of the complex crRNA and miRNA is obtained through RNAComposer. The input is an artificial sequence with the miRNA and the crRNA with secondary structure same as the complex. After fixing the complex structure it is ready for docking.

The selected server was the HDOCK server. The model that was selected after docking was the best model, meaning the model with the minimum energy. In Figure 12 is demonstrated one of the docking models. The complex crRNA/miRNA was the standard crRNA_3p with the target miR 17-3p

After evaluating the interactions of every docking model using the interactions predicted from the HDOCK server and the BIOVIA software. In Figures 13 and 14 are demonstrated the results of the interactions between the standard, mismatch, and hairpin loop complexes with the protein respectively for miR 17-3p and miR 17-5p. The interactions were studied with respect to each domain of the protein LbuCas13a and with respect to each part of the RNA Complex (Repeat, Recognition and Trigger).

It can be understood that in case of miR 17-3p, the standard design sequence interacts with the protein with the Repeat part. For the mismatch and hairpin loop designs, the interactions are formed between the Detect Part and Trigger Part of the complexes. The mismatch complex interacts solely with the NUC Lobe of the protein.

In case of miR 17-5p, the standard design complex interacts solely with REC Domain and most interactions are formed on the Repeat Part. The mismatch design complex interacts the most with REC Domain and forms only 7 interactions with the NUC Domain. Finally, the hairpin loop design interacts solely with the NUC Domain, and most of the interactions are on the Detect and Trigger part of the complex.

Standard

Mismatch

Hairpin Loop

Standard

Mismatch

Hairpin Loop

Protein Stability / crRNA - Protein Interactions Molecular Dynamics

The final step of the analysis for the CRISPR/Cas13 in terms of the crRNAs and protein system is to check the stability of the protein structure. The stability of the structure is necessary for the docking procedure to be more accurate. If the structure of a protein is unstable, then the docking algorithms will end up predicting energetically unfavorable conformations (Sponer et al., 2018).

The Molecular Dynamics simulation implemented through GROMACS software. The input file for the procedure was the fixed pdb file without the missing residues and the missing atoms.

The first step of the process was to clean the protein structure from the water. We proceeded with the creation of the topology file of the protein using the pdb2gmx module. The next steps of the process were from the basic methodology of MD simulations. The steps are the solvation, the ionization and the energy minimization of the protein. Then the next steps are the equilibration for temperature and pressure after. These steps are necessary since if they are neglected can lead to unwanted and unreliable results. Upon completion of these steps, the molecular dynamic of the protein began. The time steps chosen were 0.02ps for a total duration of 1 ns.

After running the simulation, we proceed with the calculation of the Root Mean Square Deviation of the backbone of the protein LbuCas13a, the Root Mean Square Fluctuation, Calculation of the Number of the Hydrogen Bonds formed in the protein and finally the calculation of the Solvent Accessible Surface Area of the protein. In Figures 15,16,17 and 18 are shown the plots of RMSD, RMSF, Hydrogen Bonds and SASA analysis respectively.

The Root Mean Square Deviation is representing the difference between two structures: a target structure and a reference. For the MD simulations for checking the protein stability, we evaluate the difference between the stable structure and the initial structure. In Figure 15, we can understand that the line starts to flatten even after only 400 ps. The low values of RMSD suggests that the initial structure of the protein is stable and can be used for docking purposes.

The Root Mean Square Fluctuation is a calculation of individual residue flexibility, or how much a particular residue moves (fluctuates) during a simulation. In Figure 16, we can find the residue number that demonstrates the maximum flexibility. The flexibility of a residue is revealing how feasible it is for the residue to form a bond. The more flexible a residue is, the more energetically favored is to form a stable bond. We can see that in comparison with the 2 different Lobes of the LbuCas13a protein, the REC Lobe is the more flexible hence more unstable. This agrees with the literature since the guide RNA molecule binds to the inside of the REC Lobe.

The next step of the analysis is the calculation of the Hydrogen Bonds formed inside the protein. The Hydrogen bonds are some of the most stable bonds that can be formed between 2 molecules. This is depicting their importance in biological molecules. In Figure 17, we can see that the number of hydrogen bonds after 500 ps don't not change, just fluctuates around 770. This suggests that the protein structure remains stable, since different structure results in different number of hydrogen bonds.

The final step to check the protein stability is the calculation of the Solvent Accessible Surface Area of the protein. In Figure 18, we can conclude even after 80 ps the area does not have, only fluctuates around 620 nm2. It can be easily understood that since the area does not change over time, then the surface of the structure remains the same. So, we can confirm the stability of the protein structure.

The next step for the analysis of our system is the molecular dynamics of the complex crRNA and LbuCas13a protein. For this step, we evaluate 2 different generations of crRNA for the target of miR 17-3p. The two generations of sequences are the standard design and the extra hairpin loop. For the initial coordinates of the system, we used the PDBs of the best docking models. We separated the RNA sequences from the protein and then after using the Amber force field, we used the pdb2gmx function of the GROMACS for the topology of the molecules. After the steps of the ionization and solvation, we equilibrated the system with respect to temperature and pressure. After these steps, the simulation began. The simulations of the systems were for a total time of 1 ns with time step of 1 ps.

In Figure 19 are demonstrated the plot of Root Mean Square Deviation over time for the crRNA molecules and protein. It can be concluded that in case of the interaction of LbuCas13a with the standard crRNA, the molecules are more stable since the values of RMSD are lower. On the other hand, it can be discerned that every molecule achieves a stable conformation over time.

In Figure 20 are demonstrated the Root Mean Square Fluctuation for the crRNA sequences and the LbuCas13a. Lower values of the fluctuations means more stable conformation of the specific nucleotide or residue. In the case of the protein for every simulation the residues that have lower values of RMSF are the residues that belong to the NUC lobe. For the standard design the nucleotides that belong to the repeat part exhibit the lower values of RMSF. On the other hand, the extra nucleotides of hairpin design exhibit comparable values to the repeat part, suggesting a weaker binding.

Figure 21 shows the plots of the hydrogen bonds formed in: A) LbuCas13a, B) crRNA sequences, C) between the LbuCas13a and crRNAs. The number of hydrogen bonds remain over time relatively steady suggesting that the initial conformations do not dramatically change. The initial conformations were acquired through the process of utilizing the HDOCK docking server. Overall, the number of hydrogen bonds formed between the molecules are higher for the hairpin loop design. From the calculation of the binding energy for every system, it can be shown that the binding energy in case of standard design is higher. The mean energies were calculated:

The main contribution to the energy were the Coulomb energies suggesting that the binding of the molecules was mainly due to electrostatic interactions. In Figure 22 are demonstrated the interaction energy plot for every system.

Droplet Microfluidics in-silico modelling

The field of microfluidics has recently attracted increased attention owing to its wide application in various fields of science and engineering, including biochemistry, chemistry, physics, and biotechnology. This multidisciplinary field refers to the science and technology of manipulating and allowing precise control of fluids in ultrasmall volumes (microliter, μL to femtoliter, fL), which are constrained in micrometer-scale channels.1Droplet microfluidics is one of the most important subcategories of microfluidics, which creates and manipulates discrete droplets through immiscible multiphase flows inside microchannels. Droplet microfluidics, as opposed to continuous flow microfluidics, is performed in discrete fluid droplets and this is realized by generating two-phase flows of immiscible phases, such as aqueous droplets flowing in oil. Droplet microfluidics enables the manipulation of small sample and reagent volumes in physically discrete droplets. This in turn allows high-throughput analysis without losing sensitivity.5 Droplet microfluidics samples are compartmentalized into thousands and millions of sub-samples within an immiscible phase, also known as the continuous phase (Sackmann et al., 2014, Ronshin et al., 2022).

Our research project is based on droplet microfluidic technology since it has proven that it can increase the sensitivity of the diagnostic technique and increase in general the Level of Detection (LOD). Increasing the Level of Detection constitutes feasible to detect small quantities of the biomarkers and discerning the overall overexpression. This increase is justified through the phenomenon of natural confinement effect. In general, decreasing the volume that the reactions take place, increases the concentration of the reactants and the activity of the protein. Modelling the behavior of the droplet generation is important since it can increase the efficiency of the protocol for the device manufacturing and testing (Femmer et al., 2015).

Microfluidic chip: Flow focusing Geometry 3D design

The 3D-printing technique that was selected for the device manufacturing was Digital Light Projection. Since the resolution of this technique it's limited, the designed devices should exhibit appropriate fluidic behavior without needing the channels to be less than 400 μm. With this resolution limitation the selected design was the flow-focusing devices. Generally, flow-focusing droplet generators can deliver a broader range of droplet diameters and generation rates, in comparison to the other methods of droplet formation (Van Der Linden et al., 2020).

The selected design for the microfluidic chip was a simple flow focusing geometry, with 2 oil inlets (continuous phase), 1 water inlet (disperse phase) and 1 outlet. The channel depth for all the channels selected to be 500 μm. The width of the oil inlets was 400 μm, the water inlet 500μm and the outlet channel width 700 μm. In the cross-section of the channels 4 different components that connected the water inlet with the outlet. These components had different lengths to check different types of devices. The software used for the design of these devices was SOLIDWORKS. In Figure 23 is demonstrated the channels and the overall design of the microfluidic device.

Computational Fluid Dynamics: Simulation Setup Parameters

Flow conditions play a major role in the observed performance of a flow-focusing droplet generator. The process of droplet formation is dominated by two opposing forces, the viscous forces exerted by the continuous phase, and the interfacial surface tension of the dispersed phase resisting droplet break-up. Furthermore, at low flow rates (generation rates) the effect of viscous forces becomes less dominant, while interfacial surface tension and Laplace pressure become the governing forces in dictating the process of droplet formation (Lashkaripour et al., 2021).

For the simulation of the droplet generation behavior of these devices, the software that was selected was the COMSOL Multiphysics. Different sets of inflow velocities were selected to check the flow conditions that needed to be met for the droplet formation.

The model can be shown here. Since the three dimensions constitute the simulations slow, we used the 2D model for the fluidic behavior. The physics selected for the model was the Two-Phase Flow with the Level Set method and the flow was laminar.

The geometry that was used, was only the channels with oil inlet width 400 μm and water inlet 500 μm. The outlet width was 700 μm. The length of the channels is designed to be 1 mm for the water inlet and the outlet. On the other hand, the length of the oil inlet is chosen to be 0.5 mm. 4 different geometries designed based on the middle fluidic component that connects the water inlet channel with the outlet. The length of these middle components is chosen to be 500,1000,1500 and 2000 μm.

The materials used for the simulation were water and oil. For the water material, the density chosen to be 1000 kg/m3 and the dynamic viscosity 1.0016*10-3 Pa*s. For the oil material, the density was 1850 kg/m3 and the dynamic viscosity 3.4*10-3 Pa*s. The simulation included surface tension force with coefficient 2.87*10-3 N/m.

The channel walls are chosen to be wetted with Factor of minimum element length fh 0.1. The value for the contact angle was 100 degrees. The reinitialization parameter for the level set method is set to be 0.2. The option for the mesh was finer.

The time step of the simulation is 0.001 second with a total duration of 0.8 seconds. This time duration proved to be efficient for the droplet generation.

The range of the input values for the inflow velocities are chosen to be from 0.01 m/s to 0.1 m/s for the oil and from 0.01 m/s to 0.05 m/s for every device.

Droplet Diameter and Generation Rate evaluation

The evaluation of the simulation implemented with a custom image analysis MATLAB script. The evaluation contained 2 main parts, the droplet diameter, and the Generation Rate for every flow condition. After simulation, the results were saved as picture frames of the fluidic behavior video produced by the COMSOL. After the necessary transformation, the images were converted to binaries and then the counting of the number of the droplets with its diameters began.

In Figure 25 are shown the 4 Droplet Diameter Plots with the Inflow velocities ratio, oil inflow velocity against water inflow velocity. Each plot is referring to a different geometry. Each geometry is named after the length of the middle component, 500, 1000, 1500, 2000 μm. It can be easily understood that for low inflow velocities ratio, we don't have segmented flow and no droplets are generated. We can see for the same values of inflow velocities ratio; different diameters can be achieved. The slope for every plot seems to be negative and that means that as the inflow velocities ratio is increasing, so the droplet diameter is decreasing.

500μm

1000μm

1500μm

2000μm

In Figure 26 are shown the 4 Generation Rate Plots with the Inflow velocities ratio, oil inflow velocity against water inflow velocity. Each plot is referring to a different geometry. It can be understood that when the fluidic conditions do not meet the requirements for segmented flow, the generation rate of droplets is zero. Same as the droplet diameter plots for low inflow velocities ratio values, the generation rate is zero. Overall, it seems that some specific conditions for every geometry can lead to high values of generation rate (> 100 Hz) but in most cases the generation rates are not affected by the different inflow velocities ratio.

500μm

1000μm

1500μm

2000μm

In Figures 27 and 28 are demonstrated the plots for a specific geometry, 1000 μm, for the Droplet Diameter and Generation Rate with the Inflow velocities ratio. The plots are referred to a standard value for the water inflow velocity, 0.01 -0.04. It can be understood that for the Droplet Diameter plots, the diameters of the generated droplet are clearly decreasing after increasing the inflow velocity of the oil. On the other hand, the generation rate seems to increase for the water inflow velocity condition 0.01 and 0.02, but the other conditions do not exhibit high generation rates.

From the data acquired after every simulation, the best fluidic conditions can be shown in Table 10. For every geometry 2 flow conditions were selected. The first condition is the one with the lowest droplet diameter possible. The second condition is the one with the highest generation rate possible.

Reaction Kinetics

The reactions kinetics refer to a type of mathematical modelling that through ordinary differential equations aims to predict the concentration of the reactants and the products. Using the mathematical software MATLAB, we constructed the chemical reactions of our system, and we translated them to ODEs. The main objective of this type of mathematical modelling is to predict the production of fluorescence based on the reactions' constants and the protein activity. The next step of the modelling process was to relate the experimental data with the predicted data.

In-silico modelling: Michaelis - Menten Equation

The modelling system consists of three main parts. The first part is the complex of protein LbuCas13a with the guide RNA, crRNA. The second part is the trigger RNA - miRNA target and the final part is the reporter RNA. Assuming that the protein is activated only when the guide RNA detects the target RNA, then the process can be simplified as follows.

The \( \text{Active}_{\text{Complex}} \) is referred to as the complex of the LbuCas13a/crRNA and miRNA. The reporter is the reporter RNA molecule that when the cleavage occurs, fluorescence is produced. The Intermediate state is the state of the recognition of single stranded RNA of the Active Complex and the process of the cleavage is resulting in Fluorescence and Active Complex. Based on these reactions, it is possible to study the increase of fluorescence using the Michaelis-Menten equation.

Where \( \upsilon_0 \) is the reaction velocity and \( [S] \) is the concentration of the substrate. For our system the term \( \upsilon_0 \) is the increase of the fluorescence and the substrate is the concentration of the reporter RNA. We can write our equation as:

The \( k_{\text{cat}} \) represents the rate of the enzyme activity and in our case the cleavage activity of LbuCas13a. The \( [E] \) represents the total enzyme concentration and in our case is the active protein which is the complex LbuCas13a/crRNA/miRNA.

So, our final equation will be:

The equation (4) can be modified once more if it considered the photodegradation effect with the subtraction of one more term:

Based on the equation (4) and (5), 2 basic models were constructed. The first model was with the photodegradation effect, and the figures are demonstrated in Figure 29. The time of the simulation was 10 minutes. The initial values of complex and reporter were 20nM and 1250nM respectively. For the first figure the \( r_{\text{degr}} \) chosen to be 0.01 and the values of the \( k_{\text{cleavage}} \) selected to be 1-5. For the second figure the values selected for the \( r_{\text{degr}} \) were 0.01-0.05 and \( k_{\text{cleavage}} \) to be 5. The value of \( K_m \) selected to be 500nM.

For both figures can be seen that the fluorescence is rising to a maximum that differs from every curve. For the figure A) of Figure 29, it can be discerned that every curve reaches its maximum in different time and in general lower the \( k_{\text{cleavage}} \) the peak is reached in a longer time duration. The curve with \( k_{\text{cleavage}} \) equal to 1, can be seen that when it reaches its maximum stays stable. For all the other curves, the fluorescence is decreasing. It seems that higher the \( k_{\text{cleavage}} \) the lower the final value of fluorescence will be.

For the figure B) we can discern that higher the degradation constant sooner the curve will reach its maximum. Also, it can be easily understood that the higher the degradation constant the lower the peak that will be reached for every curve. It can be shown that the fluorescence for every curve result to be 0 for \( k_{\text{cleavage}} \) 5. Finally, it is shown that the \( r_{\text{deg}} \) 0.05 has a result of a curve that decreases with the time with some fluctuations.

A

B

In Figure 30 are demonstrated the fluorescence curve with time without considering the photodegradation effect. The equation used for these curves is equation (4). In Figure 25, the parameter that changes is the cleavage constant and has values 1 to 5. We see that for the curves with 3-5 cleavage constant, they reach a maximum value of fluorescence. On the other hand, the curves with cleavage constant 1 and 2 do not reach a maximum value and the fluorescence increases steadily with time. It can be easily understood that those curves would reach the maximum value of the other if the time of the simulation was longer. The peak that is reached for every curve is based on the concentration of the reporter. This means that we will be reaching, all the reporters will be cleavaged. In Figure 31 and 32 are demonstrated the curves with cleavage rate 5 and differ the initial conditions. For Figure 31 the initial value that is changed is the trigger or miRNA concentration. This concentration means that the complex concentration will be different. We can see that every curve reaches its maximum in the same period of time. In Figure 32 are demonstrated the curves with different initial conditions of reporter RNA concentration. This means that for every curve the concentration of the substrate will differ.

The Figure 33 represents the maximum fluorescence after 10 mins for every concentration of trigger/miRNA. We can see that the fluorescence depends on concentration of the miRNA with linearity on a logarithmic scale. This linear behavior is demonstrating that different amounts of fluorescence detected can be the result of different initial values of miRNA. It is shown that the concentrations that can be detected with this technique range from 10-15 to 10-9 with small errors.

Experimental validation

The last step of the modelling process is to validate the model through the experimental results. In the next figures are demonstrated the results of the measured fluorescence with time.

1.00 μM

0.75 μM

0.50 μM

0.25 μM

In Figure 34 are demonstrated the fluorescence curves of the LbuCas13a experiments. The concentrations of the crRNA 17-3p standard, LbuCas13a and miR 17-3p selected to be 10, 20 and 1 nM respectively for all the experiments. The concentration of reporter RNA selected to have values 1, 0.75, 0.5, 0.25, 0.1 μM. In every figure it can be seen that the fluorescence values for the experiments with the miR target are higher than the control values.

Every curve shows a maximum around 10 minutes. The curves can be split into 2 main parts. The first part of the curves is the linear increasing part that the time ranges from 0 until the curves hits their maximum. The linear part is the crucial one since it shows the time needed for the detection of the target RNA. The linear behavior can be described by the mathematical formula of Michaelis - Menten enzyme kinetics which is:

Where \( dF/dt \) is the slope term of the increasing part of the curve, \( V_m \) is the maximum velocity of the reaction and \( K_m \) is the concentration of the substrate (reporter) that will result to half \( V_m \) velocity. This mathematical formula can be transformed as:

So, the first step of the analysis is to find the slopes of the linear-increasing part of the curves. Finding the maximum value of the fluorescence on each of the curves, we take all the points before the maximum and we apply the least square line fitting.

Calculating the slope of the least line square fitting, we can find the velocity for every reaction. The next step is the plot of the variables \( 1 / \upsilon_0 \) as \( Y \) and \( 1 / [\text{reporter}] \) as \( X \). Apply the least squares line fitting on the \( X \) and \( Y \), we can derive the parameters of \( K_m / V_m \) and \( 1 / V_m \) from formula (6).

After the calculation of the \( V_m \) and \( K_m \) parameters, we can calculate the catalytic activity of the enzyme based on the equation:

The \( V_m \) term is the maximum velocity of the reaction and \( E_t \) is the concentration of the active enzyme and in our case is the concentration of the target miR (\( 10^{-9} \) M).

By modifying the fluorescence Arbitrary Units to concentration, the derived values of the parameters found to be:

The second part of the curves are decreasing and hence we can apply the mathematical formula of (5):

Considering the concentration of the reporter constant this equation can be written as:

By solving this ordinary differential equation, we have the result of:

Resulting to this type of mathematical expression for the concentration of fluorescein, we can utilize an exponential fitting of this type:

After applying this fitting to the fluorescence-time data, we derived the calculated parameters and we found out that the value of \( r_{deg} \) is: