Designing optimal circuit

Logic gate circuit

The miRADAR diagnostic test uses microRNA (miRNA) biomarkers to diagnose Multiple Sclerosis (MS). However, there is no single miRNA biomarker that has the capability to detect MS, while also ruling out other diseases which might cause a change in miRNA biomarker concentration. The eight specific miRNAs that are needed to distinguish Relapsing-Remitting MS (RRMS) patients from healthy people have been found by our team in another part of the project. To effectively implement all eight of these miRNAs in our diagnostic test, toehold switches can be used to create a logic gate design, as explained in our design page. The two logic gate types that will be used in this model are AND gates and OR gates (Figure 1). We will not incorporate NOT gates for simplicity. This means we are only considering seven of the eight miRNAs found that are needed to distinguish a healthy person from an RRMS patient, since a NOT gate would be needed to incorporate the last miRNA in the circuit. For an AND gate to produce output, multiple miRNAs all need to be present at the same time, whereas an OR gate only requires one of the possible miRNAs to activate the gate.

An almost endless amount of logic gate circuit combinations are possible. However, based on our desired circuit criteria (circuit accuracy and circuit simplicity, see below), we want to find one of the most optimal circuits that can be used for a diagnostic test with seven miRNAs. We hypothesise that a circuit with more AND gates will give the highest accuracy when all miRNAs are crucial to diagnose RRMS. On the other hand, we think a circuit utilising OR gates will perform best when only a few of the found miRNAs need to be present to accurately distinguish RRMS patients from healthy persons. Experimentally, testing all possible circuits would take an incredible amount of time. However, using a computer model, many different combinations can be tested in a short time frame. Our goal was to create an algorithm which can design and test many different logic gate circuits to find the best circuit for our project. This algorithm can be used by other iGEM teams in the future if they want to create similar multi-input logic gate circuits.

Boolean Modelling

In our project, we used Boolean networks to design and test the logic gate circuit designs. Boolean networks are simple models made up of nodes and edges. Nodes are the components of a system, in our case these are the miRNAs in the sample, the toehold switch logic gates and the output that is produced. Nodes can be either in an ON or OFF state. Edges represent the interactions between different nodes of the network. These interactions function based on Boolean logic such as AND and OR relations.¹

Boolean networks offer many advantages for designing an optimal logic gate circuit. First, the logic gates themselves are already based on Boolean logic, so implementing them into a Boolean network is straightforward. Second, a Boolean network does not contain any unknown parameters, which would need to be estimated from experimental data that are absent at the design stage. Furthermore, because our diagnostic test will include a threshold detection system we assume that the miRNA concentration flowing into the logic gate circuit can only exist in two states, above the threshold (ON state) or below the threshold (OFF state). This assumption allows us to use Boolean networks for simplifying the design of our logic gate circuit.

Using a so-called update function, a Boolean network can be simulated over time. To do this, we need to provide the network with a certain initial condition in which only some parts of the network are in the ON state. The update function then uses Boolean logic, like AND and OR, to determine which nodes in the network will be in the ON/OFF state at the next iteration of the simulation. This will cause the system to end up in either of two states; one where an output is produced, or one where no output is produced. For example, in the simple logic gate circuit (Figure 1), if only miRNA 1 is present, the OR gate will go to the ON state in the next iteration, while the AND gate will stay in the OFF state and no output will be produced. To simulate a Boolean network we used an asynchronous update function. This function iterates over every order in which the nodes can updated, which ensures each possible end state is reached. The Boolean networks were simulated in Python using the PyBoolNet package.²

Circuit Optimisation

To find the best possible logic gate circuit for our diagnostic test, manually designing and testing all possible circuits would take a considerable amount of time and computational power. That’s why we created an optimisation algorithm in Python that searches for the best circuit design.

Adjacency Matrices

To build the different Boolean networks we make use of so-called adjacency matrices (Figure 2). An adjacency matrix indicates which nodes interact with each other in the Boolean network. The adjacency matrix of a circuit is always a square matrix where the number of rows and columns equal the number of nodes in the network (Figure 2). A 1 in the matrix means the node in the row is used as input in the node in the column, whereas a 0 means there is no interaction. For example, if a 1 is present at coordinate (1,4), then node 1 (row) is an input into node 4 (column). Using adjacency matrices, all interactions between logic gates and miRNAs can be easily programmed, which allows the optimisation algorithm to search through a large amount of possible circuit designs quickly.

The number of possible networks grows rapidly with number of nodes and edges. Our simple example circuit contains 6 nodes (3 miRNAs, 2 gates and 1 output, (Figure 1). This leads to 2^{36} possible unique matrices, i.e. over 60 billion networks. However, not all these networks are realistic or relevant. By setting extra criteria on our networks, most of the numbers in an adjacency matrix when designing logic gate circuits will always be 0. We can assume (i) miRNAs will never interact with each other, (ii) gates will always interact with each other in sequence (i.e., gate 2 cannot be an input into gate 1, but gate 1 can be an input into gate 2) and (iii) the last gate will always be the one producing output. These assumptions lead to only a small fraction of spaces in the adjacency matrices to be ’mutable’, i.e., the spaces where the number can be either 1 (interaction) or 0 (no interaction) (Figure 3). Only changing the mutable numbers of a circuit increases the feasibility to build a circuit optimisation algorithm, because the space in which the algorithm can search for circuits becomes a lot smaller. Following these selection criteria, our 6-node circuit now has only 2^7 = 128 possible unique matrices left to search through (Figure 3).

Testing and scoring circuits

To be able to computationally design an optimal network, the algorithm needs to be able to give a score to each circuit it designs, so that it can find the best one. In this project we created two different metrics which a circuit must oblige to: accuracy and simplicity. First, the circuit must be able to accurately predict if a specific miRNA sample should lead to MS diagnosis. To calculate the accuracy of a circuit we first determine which of the miRNAs are found to be important for diagnosing MS. We then divide these miRNAs into those considered primary miRNAs (these miRNAs must be present for an MS diagnosis) or secondary miRNAs (one of these miRNAs must be present for an MS diagnosis). Afterwards, the circuit is tested using a sample with a certain combination of miRNA. A circuit design should only produce output when all primary and one of the secondary miRNAs are present. The circuit is tested using every possible combination of miRNAs. From these results we can calculate the number of false positives (FP) and false negatives (FN) the circuit produced. The inaccuracy of the circuit design can then be calculated using the formula (FP + FN)/N_i, where N_i is the number of possible initial conditions. We need to calculate the inaccuracy instead of the accuracy because the optimisation algorithm tries to minimise the score, which in this case would be the least inaccurate circuit.

Second, the circuit must be as simple as possible to make it easier to physically build in the lab. To calculate the simplicity of the system, the total number of edges in the network is summed. Just as with the accuracy score, the algorithm will try to minimise the simplicity, which means a higher score correlates to a more complex circuit. Afterwards, a total score is calculated from the weighted sum of the complexity and the inaccuracy of the circuit. Weights of 0.1 and 0.9 are used, respectively. These weights ensure that accuracy will be the most important metric to score circuits by. In this way, the complexity of the network is only strongly considered when two circuits of differing complexity have similar accuracy.

Minimisation of circuit score

After the scoring function had been created, the PyPESTO optimisation tool was used to try and find an optimal circuit by minimising the score.³ First, a certain amount of random initial circuits was created using Latin Hypercube sampling. Latin Hypercube sampling divides the search space up into as many intervals as there are sample points and then takes a single sample from each of these intervals.⁴ This ensures the initially generated networks are relatively equally distributed across the possible interactions between nodes. Next, these circuits were optimised by flipping the 0s and 1s in the mutable part of the adjacency matrix, creating new circuit designs with higher accuracy and lower complexity. However, optimisation algorithms are made to optimise continuous values, and we are using only the binary values of 0 and 1 in our matrix. To solve this issue, the lower- and upper bounds of the values to optimise were set to -10^{-11} and 10^{-11} respectively. Afterwards, every negative value was converted into a 0 in the adjacency matrix and every positive value into a 1. This is not the only solution to the problem. If the problem would be coded in a different way, we could use other Boolean optimisation methods.⁵ However, having such a small range for the lower- and upper bound of the values to optimise ensured that the optimisation algorithm treats these values as almost binary. Every circuit designed using this technique was tested and scored afterwards. This allowed the algorithm to search for the lowest scoring circuit, which is the most optimal circuit it could find.

To run the logic gate circuit design optimisation algorithm, four parameters must be given as input. These four parameters are: the number of miRNAs important for diagnosing MS, the number of gates in the circuit design, the number of miRNAs that are considered primary (all must be present in the sample) and the number of miRNAs that are considered secondary (one of these must be present in the sample). Additionally, the number of starting points (and, as such, the number of independent optimisation runs) can be specified. A larger number of starting points increases the quality of the optimisation by lowering the average score of the optimal circuit the optimisation algorithm can find (Figure 4a-c). This decrease in score is large when using a low amount of starting points, but flattens when increasing the amount of n\_starts from an already large amount of starts. A larger number of starting points also linearly increases the computational time required to run the optimisation on a single machine (Figure 4d). This increase in computational time is more noticeable in when using a bigger number of miRNAs (Figure 4d). A consideration must be made if the gain of a lower scoring circuit is worth the extra computational time required. For our results, we used an n\_starts value of 100. This number was chosen to ensure we minimised the score following the observed trend (Figure 4), while keeping the increase in computational time manageable.

The optimisation algorithm is not able to determine the optimal number of gates that should be present in the logic gate circuit design automatically. To mitigate this, the optimisation algorithm was run multiple times with a differing number of gates to see which number logic gates would be optimal. Additionally, the algorithm cannot incorporate NOT gates into the circuit design, which limits us to using only seven of the eight miRNAs found to distinguish an RRMS patient from a healthy person. Both these limitations of the algorithm can be improved by iGEM teams in the future if they wish to design similar multi-input logic gate circuits.

Optimisation of circuit for diagnosing RRMS

As mentioned in the beginning of this page, we found seven miRNAs to be important for correctly distinguishing between an RRMS patient and a healthy person. Of these seven miRNAs, three were considered primary miRNAs and four were considered secondary miRNAs based on the importance of each miRNA in diagnosing RRMS. The optimisation algorithm explained above was run for systems with one to six logic gates. The number of starting points was set at 100 for every run.

From the results of the optimisation (Table 1) we observed that increasing the number of logic gates in the circuit design does not cause an increase in the accuracy of the circuit. However, none of the circuits obtained 100\% accuracy, which means none of the circuits was able to correctly diagnose MS in all the test cases. The best circuit the algorithm was able to find was a circuit with only a single AND gate (Figure 5).

Since miRNAs 1, 2, and 3 were considered the most important (primary) miRNAs, putting all three of these together in an AND gate can be considered the most optimal possible circuit using only one gate. The circuit should be able to make a more accurate circuit with 2 gates. This can be achieved by adding an OR gate to the circuit with the secondary miRNAs, miRNAs 4, 5, 6, and 7, as input. However, the optimisation algorithm was not able to find this circuit and instead found a sub optimal circuit consisting of a second gate which only has the first gate as input (Figure 6).

In the circuit with two logic gates (Figure 6), the second gate has no use, because it is only using the first gate and none of the miRNAs as input. Two different reasons can be given for this outcome. The first reason why only one gate is effectively used, could be because the slight increase in accuracy a second OR gate would give, is not worth the extra complexity this adds to the system according to the scoring function. Upon further investigation this was not the case. The optimal circuit found by the optimisation algorithm obtained a score of 2.1, whereas a score of 1.6 could be obtained using a 100% accurate circuit with an extra OR gate that utilises all secondary miRNAs. The 100% accurate system does, however, have a complexity value of 16 compared to the 12 found in the optimal circuit, even though the accuracy only increases from 99.21\% to 100\% (Table 1). This leads to a design choice, should accuracy be more important compared to a simple circuit with less accuracy? To change the importance of the two criteria, the weights for inaccuracy and complexity in the scoring function should be altered.

Second, it could be possible that the optimisation algorithm was not able to find a better circuit due to the number of starts being 100. Since the potential score that a 100% accurate would obtain (1.6) is lower than the minimal score the optimisation algorithm found (2.1), this was indeed true. Increasing the number of starts could lead to the lower scoring circuit being found, at the cost of increasing the computational time. Nevertheless, the logic gate circuit optimisation algorithm was able to find a simple circuit with an accuracy of 99.21\%. This circuit design can be tested in the future to find out if it is a suitable design for our miRADAR diagnostic test. These results also show that the logic gate circuit optimisation algorithm can be used and built upon by iGEM or other research teams when a design for a multi-input logic gate circuit is needed. The algorithm, including instructions on how to use it, can be found in our GitLab repository.

Dynamic modelling of toehold switch and AND gate

Ordinary Differential Equations

The logic gate circuit design algorithm finds a good circuit to use in our diagnostic test. However, the algorithm only uses a qualitative optimisation, which gives us no information on how the system behaves over time. To accurately predict and study the dynamics of the system, a quantitative model is needed. Using a quantitative model, we aim to understand the impact of leakiness on the system or which parameters have the biggest influence on the dynamics.

In our project, we designed parts for a toehold switch (Toehold Switch B) and an AND gate (two input toehold AND gate). To be able to simulate the dynamics of these designs in silico, we created Ordinary Differential Equations (ODEs). ODEs model the change in concentration of a certain system component over time. When multiple ODEs are combined, the full concentration dynamics of a system can be simulated. The ODEs of the toehold switch (model 1) and the AND gate (model 2) are defined as. \begin{align} \frac{dmR}{dt} &= 0 \tag{1a}\\ \frac{dC}{dt} &= k_b + k_c\frac{mR}{K_{M} + mR} - k_dC \tag{1b} \\ \frac{dS}{dt} &= -k_oC\frac{S}{K_{Mp} + S} \tag{1c}\\ \frac{dP}{dt} &= k_oC\frac{S}{K_{Mp} + S} \tag{1d} \end{align} \begin{align} \frac{dmR1}{dt} &= 0 \tag{2a}\\ \frac{dmR2}{dt} &= 0 \tag{2b}\\ \frac{dC}{dt} &= k_b + k_c\frac{mR1}{K_{M1} + mR1}\frac{mR2}{K_{M2} + mR2} - k_dC \tag{2c}\\ \frac{dS}{dt} &= -k_oC\frac{S}{K_{Mp} + S} \tag{2d}\\ \frac{dP}{dt} &= k_oC\frac{S}{K_{Mp} + S} \tag{2e} \end{align} In these models, we assume (i) that miRNAs do not degrade over time because they are very stable.⁶ Furthermore, we assume (ii) that there is a surplus of toehold switches in the sample, which means that the rate at which miRNAs and toehold switches form complexes is only dependent on the concentration of miRNA, the rate of complex formation k_c and the Michaelis constants K_{M1} and K_{M2}. Additionally, we assume (iii) there is no cooperative binding of miRNAs to the toehold switches, which means we can assume the model follows Michaelis Menten kinetics and no Hill coefficients are needed.⁷ Last, we assume (iv) that the output, chlorophenol red in this case, is stable and does not degrade in the time span during which we are measuring the data.⁸

The equation \frac{dC}{dt} gives us the rate at which the miRNAs and toehold switches, or AND gates, form a complex with one another. This complex can also fall apart, which is modelled using the k_d parameter. When a complex is formed, the toehold switch unfolds and the enzyme \beta-galactosidase starts to be produced. The parameter k_b controls the leakiness of the network by modelling the spontaneous unfolding of the toehold switches without any trigger present, which also produces \beta-galactosidase. When \beta-galactosidase is produced, it starts converting the yellow coloured substrate, chlorophenol red-\beta-D-galactopyranoside (CPRG), to the purple coloured product chlorophenol red, which we experimentally measured using absorbance. The rate at which substrate is converted into product is determined by the parameter k_p in the \frac{dS}{dt} and \frac{dP}{dt} equations. This assumes the formation of product follows Michealis Menten kinetics, which means there will be a decrease in product formation when the substrate concentration approaches the Michaelis constant K_{Mp}.

Data fitting

Experimental data

In wet-lab experiments performed by our team, we obtained absorbance data over time for one of our toehold switches and an AND gate. The ODE model we created can be fitted against this data to find values for all the parameters in the model. To fit the model against the data, the absorbance values first need to be converted into concentration values. This is done by subtracting the background absorbance and then using the Beer-Lambert law (A = \epsilon b C, with b=0.7cm and \epsilon=0.045\mu M^{-1}cm^{-1}) to convert the absorbance into \mu M concentration values. Since we have 3 replicates of each experiment, we can calculate the mean concentration and variance of these experiments (Figure 7) & (Figure 8).

In the experimental data graphs (Figure 7) & (Figure 8) we see that there is output (chlorophenol red) produced when the miRNA triggers are not present in both the toehold switch and the AND gate (orange and blue lines, respectively). This led us to believe there is a notable amount of leakiness in the system. To test this, the toehold switch and AND gate ODE models were fitted against both the data with miRNA triggers present and without triggers present at the same time. This way the correct value for the leakiness parameter could be found.

Parameter fitting

To build an accurate model, we must find parameter values for the model which match the experimental data. In the ODE model, the substrate consumption and product formation parts have exactly the same structure. This makes sense, because both experiments use the lacZ gene on the toehold switch to produce chlorophenol red. This means we can use the same parameter values for the \frac{dS}{dt} and \frac{dP}{dt} equations in both models. The toehold switches used in the experiments are different (Toehold Switch B and two input toehold AND gate), so the \frac{dC}{dt} part of both models will need different parameters.

To find a parameter set that can accurately match the model to the data, we first took a large Latin Hypercube sample (n=1000000) of sets for all the parameters values. This sample will take a random value for each parameter between the bounds we determined (Table 2). By taking a sample using Latin Hypercube sampling we made sure that the parameter values were equally distributed between the bounds⁴. Additionally, the parameter values were implemented on a logarithmic scale. This way we had the largest chance that a good parameter set can be found within our sample.

After generating a large amount of parameter value sets, we examined how well the model fits against the data for each parameter set. To do this, the model was simulated using the different parameter sets until t=200 minutes. The initial condition for the concentration of miRNA to simulate the model is mR(0) = 5 \mu M when the miRNA trigger is present in the sample and mR(0) = 0 \mu M when it is not present. Furthermore, we assume there is no miRNA-toehold switch complex and chlorophenol red product present at the beginning of the experiment (C(0)=0\mu M, P(0)=0\mu M). The initial concentration of substrate was calculated to be 730 \mu M in the experiment using the AND gate with both triggers and 850 \mu M in all other experiments.

When the simulation was finished, we extracted the values for the concentration of the chlorophenol red product over time and interpolated these for the eight time points in the data set (t=0, t=20, t=40, t=60, t=80, t=100, t=120 and t=180). After the interpolation we had eight values for the simulated concentration of product which we can compare to the data. We calculated a score for how well the model matches the data based on the weighted sum of squared errors between the data points and the simulation. Time point t=0 was not used in the calculation of the score, because the product concentration at t=0 will always be 0\mu M.

The method of simulating the model using a parameter set, interpolating the product concentration values at each experimentally measured time point and calculating the score was done for every parameter value set (n=1000000). One of these sets will have the lowest score and thus the simulation matches the data the best using that parameter set. To optimise the parameter values even more, we used the ’minimize’ function in SciPy.⁹ This minimisation algorithm will try to make small changes in the parameter values to minimise the score further. When the SciPy minimisation was finished we obtained a best fitting parameter value set (Table 3). The simulation of the model using the optimal parameter set against the experimental data for both the toehold switch and AND gate can be seen below (Figure 9) and (Figure 10)

Model predictions

From the results of the parameter fitting (Figure 9) and (Figure 10) we can conclude that the model replicates what we see experimentally. This means we can use the model to make predictions about the system. In the data, there is an increase in the concentration of chlorophenol red even when the miRNA trigger is not present. This can mean the system is leaky and the production of output is an effect of the leakiness as opposed to the toehold switches. To test the effect leakiness has on the system we inspected the leakiness parameter from the ODE models in more detail. The model was simulated three times for both samples (the Toehold Switch B and AND gate experiments); once with the leakiness parameter k_b included, once with the leakiness parameter removed and once without both leakiness and trigger present (Figure 11).

From these simulations we observe what happens when the leakiness parameter k_b is removed from the model. Here, the system is still producing chlorophenol red output, albeit in lower concentrations. This means that even though the system is leaky, the toehold switch and AND gate still have an effect on the chlorophenol red production and thus we can conclude that the toehold switch and AND gate are working. When both the miRNA triggers and leakiness parameters are removed, both simulations have a constant chlorophenol red concentration of 0\mu M. This is to be expected, since the non leaky system should not be producing output without any miRNA trigger present in the sample.

In the future, our ODE model can be used to make more predictions on the behaviour of the toehold switch systems. For example, a sensitivity analysis could be done on every parameter. Using sensitivity analysis, we can see the effect that increasing or decreasing the value of a parameter has on the system. This way, the model can be used to design new experiments by predicting what would happen when changing experimental conditions, such as the initial concentration of miRNA or the rate of output production by the toehold switches. When such a change has a positive effect on, for example, the speed of output production or the leakiness, we can replicate these conditions experimentally to improve our logic gate circuits. Additionally, the model can be used to study the exact dynamic behaviour of our toehold switches in more detail, increasing our knowledge about the diagnostic test overall. This way we can maximise the diagnostic value of the miRADAR test.