Team:TUDelft/Model
Modeling
The core of AptaVita is to distinguish between different vitamin concentrations by measuring the absorbance of the enzymatic product over time. Here, we simulated the temporal evolution of the product by considering the kinetics of AptaVita’s biomolecular network. The aim is to verify that our system can distinguish between different vitamin concentrations and analyze how the sensitivity of our test can be increased.
Introduction
AptaVita is a novel and modular aptazyme-based biosensor used to diagnose vitamin deficiencies from blood, contributing to increasing data availability on vitamin deficiencies and tackling hidden hunger. To provide high-quality data, our system should be able to distinguish between different vitamin concentrations. Here, we report on our deterministic kinetic model that represents the biomolecular network of AptaVita. We focused on vitamin B1 as a case study, which has a cutoff value for deficiency at 70 nM [1]. We tested our model on vitamin concentrations of 50 nM and 90 nM, representing a vitamin B1 deficient and non-deficient individual, respectively. Our model supports that our system can distinguish between different vitamin concentrations and thus suggests that the underlying design of AptaVita works as intended. Furthermore, we analyzed how the sensitivity of the AptaVita test can be improved so that smaller concentration differences can be detected. Results indicate that the aptazyme-specific parameters are of most influence and that lowering the DNA concentration also increases the sensitivity. The latter is easily adjustable in practice.
Modeling of AptaVita’s biomolecular network
Biomolecular network explanation
Fig. 1 shows a graphical representation of AptaVita’s biomolecular network. Our model simulates the temporal evolution of the product concentration in a cell-free system (CFS). The readout of the test is performed by measuring the absorbance of the product, which is proportional to the product concentration. The biomolecular network is divided into four modules:
- Gene expression (transcription-translation): The DNA template (DNA) is transcribed into a transcript consisting of the aptamer and ribozyme (the aptazyme) fused with a lacZ gene, also referred to as the uncleaved mRNA (umRNA). Upon self-cleavage of the aptazyme, the ribosome binding site of the transcript is exposed, resulting in the translation into a β-galactosidase monomer (E*). The transcription and translation steps occur in a CFS.
- Aptazyme self-cleavage: If no vitamin is bound to the aptazyme, it becomes unstable and undergoes spontaneous self-cleavage. The product is referred to as the cleaved mRNA (cmRNA).
- Vitamin binding: A vitamin molecule (Vit) reversibly binds to the aptazyme. Vitamin binding stabilizes the aptazyme, thereby preventing self-cleavage. Upon binding, an aptazyme-vitamin complex (umRNA · Vit) is formed. The bound vitamin can release spontaneously, therefore allowing the aptazyme to self-cleave spontaneously. The total vitamin concentration (bound and unbound) affects the total amount of transcripts being cleaved, thus affecting the total amount of enzyme translated and the rate of product formation.
- β-galactosidase maturation and catalysis: Maturation of a β-galactosidase monomer results in a functional β-galactosidase tetramer (E). The functional enzyme catalyzes the conversion of the substrate chlorophenol-β-D-galactopyranoside (CPRG) (S) into the product chlorophenol red (CPR) (P).
AptaVita’s read-out is summarized as follows: the higher the vitamin concentration in a sample, the slower the product formation (and increase of absorbance), and vice versa.
Cell-free gene expression
Gene expression is an integral module of the biomolecular network. Hence, we performed a literature study to find a mathematical model that properly describes gene expression in a CFS. Stögbauer et al. expressed the reporter protein Green Fluorescent Protein (GFP) in the PURExpress CFS and they developed a deterministic model describing the resulting expression profile of GFP [2]. We qualitatively validated their gene expression model by expressing a construct containing the reporter protein Yellow Fluorescent Protein (YFP) in the PURExpress CFS. This YFP-construct did not contain an aptazyme sequence. The measured YFP fluorescence over time is shown in Fig. 2.
The expression profile of YFP in the PURExpress CFS (Fig. 2) matches qualitatively those measured by Stögbauer et al. [2]. Three different regimes, which are characteristic for CFSs, were identified in the expression profile [3]. The first regime (i) is the transient regime, during which gene expression is initiated. The second regime (ii) is the steady-state regime, where the rate of gene expression is constant over time. The third regime (iii) is the plateau regime, where the gene expression stops due to depletion of resources [4]. These results substantiate the use of the Stögbauer et al. model for the gene expression module of our network.
Modeling gene expression in a cell-free system
For an overview of all assumptions and species involved in the model, see Tab. 1 and Tab. 2, respectively.
Building upon an existing deterministic cell-free expression model
We adapted the minimalistic rate equation model for protein expression in the PURExpress CFS by Stögbauer et al. [2] and extended it by including the aptazyme self-cleavage, vitamin binding, and β-galactosidase kinetics. The following set of differential equations describes our model:
Eq. 1, 2, 6, and 7 describe the gene expression (transcription-translation) kinetics. Eq. 2 and 7 were directly adopted from Stögbauer et al., whereas Eq. 1 and 6 were adjusted to include species of AptaVita's biomolecular network [2]. Transcription and translation are described by Michaelis-Menten like equations to account for saturation effects. The prefactors TsR and TlR account for expression saturation by representing finite transcription- and translation resources, respectively. TsR consumption is dependent on the DNA concentration and are consumed during the transcription process with scaling factor kcs. This is in contrast to the TlR, which depletes independently from the DNA concentration described by δTlR. Furthermore, kts and ktl are the rate constants of transcription and translation, respectively. Ks, Kl, and KTlR are the Michaelis constants for transcription and transcription resources, translation, and translation resources, respectively. δmRNA is the mRNA degradation rate constant for umRNA and cmRNA [2].
Eq. 5 describes the aptamer self-cleavage kinetics, where kc is the cleaving rate constant of the aptazyme. Furthermore, Eq. 3 and 4 describe the vitamin binding kinetics. Here, kon and koff are the rate constants of the vitamin binding to- and releasing from the aptamer, respectively.
Lastly, Eq. 8, 9, and 10 describe the enzyme maturation and catalysis kinetics. Here, kmat is the maturation rate constant of β-galactosidase. The factor ¼ in front of the maturation rate constant is required because β-galactosidase is a tetramer. kcat and Km are the catalytic rate constant and Michaelis constant of β-galactosidase with CPRG as substrate, respectively.
Model parameters
An overview of all parameters is shown in Tab. 3. The initial conditions used in our model were a DNA concentration of 3 nM and a S concentration of 1000 μM, unless stated otherwise. The latter concentration was based on an estimation of the maximum product concentration for which the absorbance stays in the linear range [5]. Furthermore, the initial free Vit concentrations were set at 50 and 90 nM, representing a deficient and non-deficient sample, respectively [1]. The prefactors TsR and TlR were set to 1 [2]. Stögbauer et al. optimized their model parameters by fitting measured time courses of protein expression, and have shown that their parameter fits are robust [2]. Hence, we used their values where possible or adapted them to our system. The parameters needed to extend the model were obtained elsewhere from the literature.
The kts was adjusted to the length of our transcript, which is 3297 base pairs [6]. Stögbauer et al. reported a kts of 2.2 NTP s-1 per T7 RNA polymerase. The concentration of T7 RNA polymerase in PURExpress is approximately 100 nM [2]. The adjusted kts is equal to 6.7 · 10-5 μM s-1 and was calculated as follows:
Similarly, the ktl was adjusted to the length of the β-galactosidase monomer, which is 1023 amino acids [7]. Stögbauer et al. reported a ktl of 0.03 AA s-1 per ribosome. The concentration of ribosomes in PURExpress is approximately 2.4 μM [2]. The adjusted ktl is equal to 7.0 · 10-5 μM s-1 and was calculated as follows:
Furthermore, the value for koff was based on reported values by Townsend et al. [8]. Townsend et al. measured koff values for aptazymes generated with the help of the DRIVER method, the method for our in vitro evolution experiments. The dissociation constant (KD) was determined using the Aptagen database [9]. This database contains KD values for aptamers binding small molecules. Based on these values, we chose a KD of 50 nM and calculated kon from KD and koff as follows:
Moreover, the value for kc was obtained from Kvorova et al., as this rate constant was measured for a catalytic core - the specific location within the secondary structure where cleavage occurs - which was the same as our ribozyme [10].
Lastly, the value for kmat was obtained from Nichtl et al., who identified the structural rearrangement of β-galactosidase dimers to be the rate-limiting step in tetramer formation [11]. Hence, the value for this rate constant was taken as value for kmat. The kcat was obtained from a model fit of the β-galactosidase kinetics with CPRG as the substrate and the Km was taken from a Michaelis-Menten assay with CPRG as the substrate [8, 12].
In this section, an overview of all the assumptions is shown (Tab. 1).
| Assumption | Consequence for the model |
|---|---|
| Michaelis-Menten kinetics occurs under quasi steady-state [13] | No consideration of the enzyme-substrate complex |
| No inhibitory effects of any species | No inhibition terms are included |
| Vitamin binding prevents aptazyme degradation | Degradation only occurs for unbound umRNA and cmRNA |
| No degradation of DNA or proteins [2] | No degradation term for DNA, E*, and E |
| All reaction components are present in high copy numbers [2] | A deterministic approach is appropriate |
| No evaporation | Constant reaction volume |
| No diffusion limitations | Reaction rates are rate-limiting |
In this section, an overview of all the species used in our model is shown (Tab. 2).
| Species | Description |
|---|---|
| DNA [μM] | Template DNA |
| umRNA [μM] | Uncleaved mRNA (aptazyme fused with lacZ gene) |
| TsR [-] | Transcription resources |
| Vit [μM] | Total free vitamins (B1) |
| umRNA · Vit [μM] | Uncleaved mRNA-vitamin complex |
| cmRNA [μM] | Cleaved mRNA (lacZ gene) |
| E* [μM] | Monomer subunit of enzyme (β-galactosidase monomer) |
| TlR [-] | Translation resources |
| E [μM] | Mature enzyme (β-galactosidase tetramer) |
| S [μM] | Substrate (CPRG) |
| P [μM] | Product (CPR) |
In this section, an overview of the parameter constants is shown (Tab. 3).
| Parameter constant | Description | Value |
|---|---|---|
| kts | Transcription rate constant | 6.7 · 10-5 [μM s-1] [2, 6] |
| Ks | Michaelis constant of transcription and TsR | 8.5 · 10-3 [μM] [2] |
| koff | Off-rate constant of Vit releasing from the aptamer | 1 · 10-2 [s-1] [9] |
| kon | On-rate constant of Vit binding to the aptamer | 2 · 10-1 [μM s-1] [10] |
| kc | Cleaving rate constant of the aptazyme | 1.7 · 10-2 [μM] [10] |
| δmRNA | Degradation rate constant of umRNA and cmRNA | 1.3 · 10-5 [s-1] [2] |
| kcs | Scaling factor for TsR | 1.8 · 10-4 [s-1] [2] |
| ktl | Translation rate constant | 7.0 · 10-5 [μM s-1] [2, 7] |
| Kl | Michaelis constant of translation | 6.58 · 10-2 [μM] [2] |
| kmat | Maturation rate constant of E* | 0.5 · 10-3 [s-1] [11] |
| δTlR | Degradation rate constant of TlR | 7.5 · 10-5 [s-1] [2] |
| KTlR | Michaelis constant of TlR | 6 · 10-5 [s-1] [2] |
| kcat | Catalysis rate constant of E | 5.14 · 101 [s-1] [8] |
| Km | Michaelis constant of S | 50 [μM] [12] |
Model implementation
We implemented our model in Python, as it is an open-source programming language and thus contributes to the values of iGEM. For our first implementation, we programmed the model in a combination of pure Python and NumPy, an often-used library for scientific computing [14]. When we benchmarked the runtime of our implementation, running our model for 7200 seconds with timesteps of 0.01 seconds took 16.9 seconds. To improve the performance, we used the Numba library. This Python library has a high-performance Python just-in-time compiler that compiles certain Python code into machine code. This allows for significant runtime improvements and allows for C, C++, and Fortran-like performance in some instances [15]. After including the compilation in our Python code, we benchmarked the same model run with the same total time and timesteps of 0.01 seconds. The model now took 0.0392 seconds, which is a 431 fold performance improvement. The use of the Numba library was crucial, as it allowed us to run the simulations often and with a diverse range of parameter values, which was needed for the sensitivity analyses we performed. A result created with our model implementation is shown in Fig. 3.
It is worth noting that Numba is not able to speed up all possible Python code. Without proper care, Numba can actually slow down an implementation. Numba performs well with numerically oriented code and code with many NumPy arrays and/or for loops. This was the case for our model implementation, where many calculations need to happen for a large number of loops. However, Numba does not perform well with Python objects such as lists or dictionaries. Therefore, we would recommend future iGEM teams to use the Numba library if their code is numerically oriented, since proper use can lead to significant runtime improvements without major adjustments to existing Python code.
To perform the sensitivity analysis we used the SALib library, which is an open-source Python library that contains implementations of often-used sensitivity analysis methods. This library was used for Morris method trajectory generation and for the analysis of the results [16]. The source code of our implementation and figures can be found here.
Sensitivity analysis
The underlying principle of the AptaVita test is that a change in vitamin concentration leads to a change in the rate at which product is formed. By showing that the product formation is sensitive to changes in the vitamin concentration, we can validate that our system works as intended. To this end, we performed a sensitivity analysis. A sensitivity analysis establishes the degree of sensitivity of the model output with respect to specific input changes, such as the vitamin concentration. The analysis also provided information on the robustness of the model. Besides validating the technical design of the AptaVita test, we analyzed how the sensitivity of the test can be improved. A second sensitivity analysis showed possible adjustments of our system, so that smaller differences in vitamin concentrations can be detected.
We used the Morris method for both sensitivity analyses [17]. We chose this method because of its capability to handle non-linearity in the model, and interactions between factors. This is important, as our system is non-linear and varying factors simultaneously may affect the output to a larger extent than varying them individually. Additionally, the Morris method performs well at a relatively low computational cost compared to other methods that can handle non-linearities and interactions, such as the Sobol method [18].
To perform the Morris method, it is required to specify a range of values for all the factors that one wants to determine the sensitivity of. The ranges and the parameters that we used for the Morris analysis can be found here. The Morris method outputs sensitivity measures, of which we used the μ and μ*. The sensitivity measure μ estimates the average impact on the output of varying a factor from the lower to the upper bound of its range. A positive impact indicates that increasing a factor will lead to an increase in the output. A negative impact indicates that increasing a factor will lead to a decrease in the output. It could happen that positive and negative impacts cancel each other out and that μ is too low to properly reflect the impact. To prevent this, μ* is also calculated, which is the average of the absolute values of the impact. An explanation of the Morris method itself can be found below.
Assume that the sensitivity of k input variables, which are initial conditions or parameters, needs to be determined. For each variable, a range of values is defined. These variables are called factors. Each range is discretized into a grid consisting of p points, i.e. levels [17]. Thus, when the sensitivity analysis is done with k factors, and each range is split up into p levels, the factor space contains k dimensions and p different levels along each axis.
For the analysis, a random starting point in the factor space is picked. The model output is calculated for the corresponding combination of factor values. One of the factor values is changed to a level that is two levels away from the starting point. This can be visualized as a step in the factor space, for example from the first level to the third level. The difference between the new output and the old output is calculated. This difference is normalized with respect to the step size, so that the impact is extrapolated over the full range of the factor. The resulting value is an elementary effect of the corresponding factor. Subsequently, the value of a different random factor is changed, and the elementary effect of this new step is calculated. This process is repeated by taking new steps until each factor has been varied once. All these different steps together form a trajectory [17]. The more trajectories are used, the better the factor space is sampled. However, this also increases the computational time and memory requirements. Therefore, it is advantageous to first generate a high number of random trajectories. From these, r optimal trajectories are selected. Only these are used in the analysis. The optimal trajectories are chosen such that the distance between the trajectories is maximized. In this way, a larger part of the factor space is sampled than by using the same amount of random trajectories. In our cases, we generated 1000 random trajectories and selected r = 70 optimal trajectories from these [19]. The number of simulations that is needed for the analysis is equal to r(k+1) [18]. A visualization of the Morris method is shown in Fig. 4.
For each factor, sensitivity measures are calculated from the elementary effects of that factor, namely μ, μ*, and the standard deviation (σ). Both μ and μ* are measures for the impact of a factor on the model output when the factor is varied from the lower to the upper bound of its full range. The Morris method is suited to differentiate between factors that have a large influence on the output and those that have not. It is not quantitative enough to rank the factors in order of most impact [17]. Information about the non-linearity and interactions between factors is contained in σ. These cannot be distinguished, as both will result in the elementary effect being different in other parts of the factor space. In the case of non-linearity, the elementary effect will differ depending on the level of the factor itself, and in the case of interactions they will differ depending on the levels of the other factors.
Validation of the underlying design of the AptaVita test
Introduction
The first sensitivity analysis was performed to determine whether the product formation is sensitive to the total vitamin concentration. Showing that this is the case suggests that our system works as intended. The sensitivity analysis also provides information about the robustness of the model and whether we can make further qualitative statements about our system. The results are shown in Fig. 5. The ranges and the parameters that were used for the first Morris analysis can be found here.
Results
The factors can be grouped in three categories, based on the sensitivity of the product formation to that factor (Fig. 5). The 95%-bootstrapped confidence intervals of the μ* of the categories do not overlap, except for the beginning and end. The first category consists of the factors to which the product formation is not sensitive, as the μ* peak values are lower than 15 μM. This is low compared to the maximal product output of 1000 μM. The second category consists of the factors to which the product formation is sensitive. The μ* peak values are between 35 and 110 μM. Since the maximal product output is 1000 μM, varying these factors will have a noticeable effect on the product formation. The third category consists of the factors to which the product formation is very sensitive. This is shown by the fact that the μ* peak values are higher than 250 μM, which is relatively high compared to the maximal product output of 1000 μM.
Conclusions
The total vitamin concentration belongs to the sensitive category. This shows that changes in vitamin concentration lead to changes in product formation, which is the basis of our AptaVita test. This result is illustrated with the vitamin B1 case study in Fig. 6, where it can be clearly seen that the product formation curve for a total vitamin concentration of 50 nM is different from the curve with a total vitamin concentration of 90 nM.
Most of the other factors belong either to the sensitive or the very sensitive category, indicating that the model output is not robust to changes in these factors. The factors belonging to the sensitive group are mostly aptazyme parameters. For the purpose of a qualitative understanding of our system, it is not essential that these parameters have values corresponding to an existing aptazyme, as long as these values could be correct for a potential aptazyme. The factors belonging to the very sensitive category are mostly parameters for gene expression. These parameter values were taken and adjusted accordingly from Stögbauer et al., who have shown their parameter fits to be robust, indicating that our adjusted parameters are likely approximately correct, allowing for further qualitative analysis [2].
Our CFS is lyophilized, which may influence parameter values. Pardee et al. have shown for a different CFS that their lyophilized CFS showed activity comparable to solution-phase reactions, the type studied by Stögbauer et al. [2, 20]. Hence, we postulate that the lyophilization of the CFS does not have a significant effect on parameter values. However, we would recommend measuring the values of the parameters in the wet lab under the conditions and with the aptazyme of the final AptaVita test to improve the model validity.
Increasing discernibility between vitamin concentrations
Introduction
Our model shows that the total vitamin concentration affects the product concentration over time. However, the difference between the curves for a total vitamin concentration of 50 nM and 90 nM, corresponding to the vitamin B1 case study, is rather small (Fig. 7a). To improve the sensitivity of our test, we investigated how the discernibility between different vitamin concentrations can be improved. To this end, we performed a second sensitivity analysis. The ranges and the parameters used for the second Morris analysis can be found here. The area between two different product curves was chosen as a metric for discernibility. The area can be physically interpreted as the difference in product concentration summed over time. The larger the area, the better two curves can be distinguished (Fig. 7b).
The resulting μ and μ* of a factor both estimate the change in the area by varying this factor from the lower to the upper bound of the range. The new area can be calculated by adding μ to the reference area. The same can be done for μ*. For the interpretation of the results, it is convenient to normalize the new area to the reference area. This gives a fold change of the area, since it describes how much the area changes between the reference model run and the subsequent model run. For example a fold change of 1.2 indicates that the area increased by 20% and a fold change of 0.8 indicates that the area decreased by 20%. The reference area had as initial conditions a S concentration of 1000 μM, a DNA concentration of 3 nM, and a Vit concentration of 50 nM and 90 nM. The results are shown in Fig. 8. The fold change is calculated as follows:
Results & Conclusion
According to Fig. 8, the area is sensitive to changes in aptazyme-associated parameters. The area is also sensitive to most of the other factors. However, most factors are either hard to adjust or not adjustable, since they are inherent to the CFS or the enzyme. An important exception is the DNA concentration, for which the fold change of the area for μ* was 1.76. The DNA concentration can be adjusted in the laboratory. Therefore, the μ of the DNA concentration suggests that we can practically modify the sensitivity of our test by adjusting the DNA concentration. However, it should also be taken into account that lowering the DNA increases the total time that is required for the test, as protein expression and thereby substrate conversion are slowed down. A test is only classified as rapid according to the REASSURED criteria of the World Health Organization if it can be completed within 15 to 60 minutes [21]. Additionally, too low concentrations (<100 pM) might cause the gene expression to become unstable amongst different samples. Therefore, we would recommend testing out different DNA concentrations for the AptaVita system and measuring the variability. By doing so, an optimum can be found between the total test time, gene expression stability, and the test sensitivity.
As mentioned above, the area is most sensitive to changes in aptazyme-associated parameters, shown by kon, koff and kc having a fold changes for μ* of 3.25, 3.56, and 3.58, respectively. As these values are the most extreme, this suggests that improving aptazyme characteristics may have the most impact on the test sensitivity. Thus, it can be of value to study these further. The area being sensitive to changes in kon, and koff, makes sense intuitively, since the kon and koff determine the dissociation constant of the aptazyme KD (= koff / kon). A lower KD will lead to a higher concentration of umRNA · Vit complexes, and binding of the vitamin to the aptazyme inhibits cleavage. This lowers the translation rate, and indirectly lowers the rate of product formation. Also lowering the kc has been shown to increase the sensitivity significantly. To gain a deeper understanding of the kc on the system and a better understanding of the influence of KD, we looked at the formation of umRNA and cmRNA over time, which is discussed here.
The values of kon, koff, and kc are dependent on the outcome of the in vitro evolution of the aptazyme. Although the aptazyme parameter values cannot be as easily adjusted as the DNA concentration, we hypothesize that by adapting the original DRIVER method, we can bias the aptazymes to increase the KD and lower the kc. By performing relatively more binding cycles compared to cleaving cycles, aptazymes could be biased towards a lower cleaving rate constant, since this will be evolutionarily advantageous. In addition, the kc may also be lowered with the help of specific point mutations in one of the aptazyme loops. Point mutations at specific positions have been shown to decrease the cleaving rate constant by a factor of 10 [10]. However, these point mutations could have a converse effect, since they might decrease the binding affinity of the aptazyme to the target ligand. To illustrate the effect that lowering the kc would have, an animation of the product concentration over time with a varying the kc is shown in Fig. 9.
To gain a better understanding of kc and KD their influence on the sensitivity of the test, and thereby of kon and koff, we looked at the formation of umRNA and cmRNA over time. The concentration of umRNA and umRNA · Vit both increase during the initial phase of the test. As the concentration of umRNA increases, the amount of umRNA that cleaves increases as well, since the cleavage is governed by a first order reaction (Eq. 5). This will continue until an approximate equilibrium is reached at which the production of umRNA is about equal to the formation of cmRNA. For clarity, we call this point the production-cleavage equilibrium. The production-cleavage equilibrium point is at the time point where the difference between the umRNA formation rate and the cmRNA formation rate is the smallest. From this point onwards the rate of change in cmRNA starts to follow that of umRNA production. Most of the difference in cmRNA concentration builds up before the production-cleavage equilibrium, as the difference in the rate of change of cmRNA is much larger before the equilibrium is established (Fig. 10). The later the production-cleavage equilibrium is established, the larger the difference in cmRNA concentration for different vitamin concentrations will be. These differences in cmRNA are amplified in the following translation and enzymatic step, leading to a higher sensitivity of the test.
Lowering the KD delays the time at which the production-cleavage equilibrium is reached. If the KD is lowered, the affinity of the aptazyme is increased and thus more of the aptazymes will be bound at the same total vitamin concentration. Thus, a decrease in KD will lead to an increase in the concentration of umRNA · Vit. As binding of the vitamin to umRNA prevents cleavage, the formation rate of cmRNA decreases. As a result, it will take longer until the rate of change of cmRNA is approximately equal to the production rate of umRNA. So the production-cleavage equilibrium is delayed. As the area is very sensitive to changes in kon en koff, and KD is the fraction of the two, this effect will be relatively large.
Lowering kc also delays the production-cleavage equilibrium. As kc is lowered, less of the umRNA will be cleaved per unit time, and thus a higher umRNA concentration is needed for the formation of cmRNA to equal the production umRNA. The model also shows this, as the peak of the umRNA concentration increases for lower values of kc (Fig. 11). Additionally, the effect of lowering kc is strengthened by vitamin binding. Since the umRNA concentration is higher, more umRNA will be bound to a vitamin. This again decreases the concentration of umRNA, and thereby the formation of cmRNA.
Discussion
On this page we reported on a deterministic model representing the biomolecular network of AptaVita. During our modeling effort, we focused on vitamin B1 as a case study and tested our model on vitamin concentrations of 50 nM and 90 nM, representing a vitamin B1 deficient and non-deficient individual, respectively. We performed multiple sensitivity analyses to make qualitative statements about our system. We showed that the AptaVita systems allows for the discrimination between vitamin concentrations such as those of the B1 case study. We have also analyzed how the sensitivity of our test can be improved. Lowering the DNA template concentration increases the sensitivity in an easily adjustable manner. Furthermore, improving aptazyme characteristics also leads to a significant improvement in the sensitivity. Either increasing kon or decreasing koff or kc improves the sensitivity of the test. We have touched upon methods to improve these characteristics. However, most of these adjustments will lead to an increase in test time. For implementation in the field, an optimum between test sensitivity and test duration would have to be found.
All information for the sensitivity analysis with the product concentration as output. Morris parameters:
- Sampled trajectories = 1000
- Optimal trajectories = 70
- Levels = 4
The model output is the product concentration at every time point. The sensitivity is measured at every time point. The parameter constants were varied as follows: parameter constant · ⅓, parameter constant · 3. For exact values see below.
| Parameter constant | Lower bound value | Upper bound value |
|---|---|---|
| kts | 2.23 · 10-5 [μM s-1] | 2.01 · 10-4 [μM s-1] |
| Ks | 2.88 · 10-3 [μM] | 2.55 · 10-2 [μM] |
| koff | 3.33 · 10-3 [s-1] | 3.00 · 10-2 [s-1] |
| kon | 6.67 · 10-2 [μM s-1] | 6.00 · 10-1 [μM s-1] |
| kc | 5.67 · 10-3 [μM] | 5.10 · 10-2 [μM] |
| δmRNA | 4.33 · 10-6 [s-1] | 3.90 · 10-5 [s-1] |
| kcs | 6.00 · 10-5 [s-1] | 5.40 · 10-4 [s-1] |
| ktl | 2.33 · 10-5 [μM s-1] | 2.10 · 10-4 [μM s-1] |
| Kl | 2.19 · 10-2 [μM] | 1.97 · 10-1 [μM] |
| kmat | 1.67 · 10-4 [s-1] | 1.50 · 10-3 [s-1] |
| δTlR | 2.50 · 10-5 [s-1] | 2.25 · 10-4 [s-1] |
| KTlR | 2.00 · 10-5 [s-1] | 1.80 · 10-4 [s-1] |
| kcat | 1.71 · 101 [s-1] | 5.14 · 102 [s-1] |
| Km | 16.7 [μM] | 150 [μM] |
| Initial condition | Lower bound value | Upper bound value |
|---|---|---|
| DNA* | 3.00 · 10-4 [μM s-1] | 6.00 · 10-3 [μM s-1] |
| S** | 1000 [μM] | 1000 [μM] |
| Vit | 2.33 · 10-2 [μM] | 2.10 · 10-1 [μM] |
*The DNA concentration was varied over a larger range, since the DNA concentration can be easily adjusted in the final system.
**The substrate concentration was not varied.
The simulation was run for 7200 seconds with timesteps of 0.01 seconds.
All information for the sensitivity analysis with the area as model output. Morris parameters:
- Sampled trajectories = 1000
- Optimal trajectories = 70
- Levels = 4
The model output is the area between product curve with a vitamin concentration of 50 nM and product curve with a vitamin concentration of 90 nM. The parameter constants were varied as follows: parameter constant · ⅓, parameter constant · 3. For exact values see below.
| Parameter constant | Lower bound value | Upper bound value |
|---|---|---|
| kts | 2.23 · 10-5 [μM s-1] | 2.01 · 10-4 [μM s-1] |
| Ks | 2.88 · 10-3 [μM] | 2.55 · 10-2 [μM] |
| koff | 3.33 · 10-3 [s-1] | 3.00 · 10-2 [s-1] |
| kon | 6.67 · 10-2 [μM s-1] | 6.00 · 10-1 [μM s-1] |
| kc | 5.67 · 10-3 [μM] | 5.10 · 10-2 [μM] |
| δmRNA | 4.33 · 10-6 [s-1] | 3.90 · 10-5 [s-1] |
| kcs | 6.00 · 10-5 [s-1] | 5.40 · 10-4 [s-1] |
| ktl | 2.33 · 10-5 [μM s-1] | 2.10 · 10-4 [μM s-1] |
| Kl | 2.19 · 10-2 [μM] | 1.97 · 10-1 [μM] |
| kmat | 1.67 · 10-4 [s-1] | 1.50 · 10-3 [s-1] |
| δTlR | 2.50 · 10-5 [s-1] | 2.25 · 10-4 [s-1] |
| KTlR | 2.00 · 10-5 [s-1] | 1.80 · 10-4 [s-1] |
| kcat | 1.71 · 101 [s-1] | 5.14 · 102 [s-1] |
| Km | 16.7 [μM] | 150 [μM] |
| Initial condition | Lower bound value | Upper bound value |
|---|---|---|
| DNA* | 3.00 · 10-4 [μM s-1] | 6.00 · 10-3 [μM s-1] |
| S** | 1000 [μM] | 1000 [μM] |
| Vit*** | - | - |
*The DNA concentration was varied over a larger range, since the DNA concentration can be easily adjusted in the final system.
**The substrate concentration was not varied.
***The vitamin concentration was not varied. The vitamin concentration is fixed at 50 nM and 90 nM.
The simulation was run for 10800 seconds with timesteps of 0.01 seconds.
References
- Albaugh, V.L., Williams, D.B., Aher, C.V., Spann, M.D., & English, W.J. (2021). Prevalence of thiamine deficiency is significant in patients undergoing primary bariatric surgery. Surgery for Obesity and Related Diseases, 17(4), 653-658. https://doi.org/10.1016/j.soard.2020.11.032
- Stögbauer, T., Windhager, L., Zimmer, R., & Rädler, J.O. (2012). Experiment and mathematical modeling of gene expression dynamics in a cell-free system. Integrative Biology, 4(5), 494-501. https://doi.org/10.1039/c2ib00102k
- Doerr, A., De Reus, E., Van Nies, P., Van der Haar, M., Wei, K., & Kattan, J. Wahl, A., & Danelon, C. (2019). Modelling cell-free RNA and protein synthesis with minimal systems. Physical Biology, 16(2), 025001. https://doi.org/10.1088/1478-3975/aaf33d
- Marshall, R., & Noireaux, V. (2019). Quantitative modeling of transcription and translation of an all-E. coli cell-free system. Scientific Reports, 9(1). https://doi.org/10.1038/s41598-019-48468-8
- Sigma-Aldrich (2018). Certificate of analysis: Chlorophenol red. Adopted from https://www.sigmaaldrich.com/certificates/COFA/19/199524/199524-BULK_______MKCG3138__.pdf on 19/10/2021.
- Macdonald, L.E., Durbin, R.K., Dunn, J.J., & McAllister, W.T. (1994). Characterization of two types of termination signal for bacteriophage T7 RNA polymerase. Journal of Molecular Biology, 238(2), 145-158. https://doi.org/10.1006/jmbi.1994.1277
- Juers, D.H., Matthews, B.W., & Huber, R.E. (2012). LacZ β‐galactosidase: Structure and function of an enzyme of historical and molecular biological importance. Protein Science, 21(12), 1792-1807. https://doi.org/10.1002/pro.2165
- McNerney, M.P., Zhang, Y., Steppe, P., Silverman, A.D., Jewett, M.C., & Styczynski, M.P. (2019). Point-of-care biomarker quantification enabled by sample-specific calibration. Science Advances, 5(9), eaax4473. https://doi.org/10.1126/sciadv.aax4473
- Aptagen. (2019). Aptamer Database. Adopted from https://www.aptagen.com on 06/10/2021.
- Khvorova, A., Lescoute, A., Westhof, E., & Jayasena, S.D. (2003). Sequence elements outside the hammerhead ribozyme catalytic core enable intracellular activity. Nature Structural & Molecular Biology, 10(9), 708-712. https://doi.org/10.1038/nsb959
- Nichtl, A., Buchner, J., Jaenicke, R., Rudolph, R., & Scheibel, T. (1998). Folding and association of β-galactosidase. Journal of Molecular Biology, 282(5), 1083-1091. https://doi.org/10.1006/jmbi.1998.2075
- Prasad, V.V., & Fliesler, S.J. (1994). Identification of β-galactosidase activity in purified bovine retinal rod outer segments. Current Eye Research, 13(5), 377-384. https://doi.org/10.3109/02713689409167302
- Kang, H.W., KhudaBukhsh, W.R., Koeppl, H., & Rempała, G.A. (2019). Quasi-steady-state approximations derived from the stochastic model of enzyme kinetics. Bulletin of Mathematical Biology, 81(5), 1303-1336. https://doi.org/10.1007/s11538-019-00574-4
- Harris, C., Millman, K., Van der Walt, S., Gommers, R., Virtanen, P., & Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N., Kern, R., Picus, M., Hoyer, S., Van Kerkwijk M., Brett, M., Haldane, A., Fernández Del Río, J., Wiebe, M., Peterson, P., Gérard-Marchant, P., Sheppard K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., & Oliphant, T. (2020). Array programming with NumPy. Nature, 585(7825), 357-362. https://doi.org/10.1038/s41586-020-2649-2
- Lam, S., Pitrou, A., & Seibert, S. (2015). Numba: a LLVM-based Python JIT compiler. LLVM 15’: Proceedings of the second workshop on the LLVM compiler infrastructure in HPC. 7(6). https://doi.org/10.1145/2833157.2833162
- Herman, J., & Usher, W. (2017). SALib: An open-source Python library for sensitivity analysis. The Journal of Open Source Software, 2(9), 97. http://doi.org/10.21105/joss.00097
- Saltelli, A., Ratto, M., & Andres, T. (2008). Global sensitivity analysis: the primer. Hoboken, N.J.: Wiley, 109-128.
- Qian, G., & Mahdi, A. (2020). Sensitivity analysis methods in the biomedical sciences. Institute of Biomedical Engineering, Oxford, U.K.
- Ruano, M., Ribes, J., Seco, A., & Ferrer, J. (2012). An improved sampling strategy based on trajectory design for application of the Morris method to systems with many input factors. Environmental Modelling & Software, 37, 103-109. https://doi.org/10.1016/j.envsoft.2012.03.008
- Pardee, K., Green, A., Ferrante, T., Cameron, D., DaleyKeyser, A., Yin, P., & Collins, J. (2014). Paper-based synthetic gene networks. Cell, 159(4), 940-954. https://doi.org/10.1016/j.cell.2014.10.004
- Land, K., Boeras, D., Chen, X., Ramsay, A., & Peeling, R. (2018). REASSURED diagnostics to inform disease control strategies, strengthen health systems and improve patient outcomes. Nature Microbiology, 4(1), 46-54. https://doi.org/10.1038/s41564-018-0295-3